diff --git a/doc/src/bond_special.rst b/doc/src/bond_special.rst index f7dc43a1b2..5f2cebbb5a 100644 --- a/doc/src/bond_special.rst +++ b/doc/src/bond_special.rst @@ -73,7 +73,7 @@ ensure that the new bonds created by this style do not create spurious Specifically 1-2 interactions must have weights of zero, 1-3 interactions must either have weights of unity or :doc:`special_bonds angle yes ` must be used, and 1-4 interactions must -have weights of unity or :doc:`special_bonds dihedral yes ` +have weights of unity or :doc:`special_bonds dihedral yes ` must be used. If this command is used to create bonded interactions between @@ -95,7 +95,7 @@ compute interactions for individual pairs of atoms. Manybody potentials are not compatible in general, but also some other pair styles are missing the required functionality and thus will cause an error. -This command is not compatible with long-range Coulombic interactions. If a +This command is not compatible with long-range Coulombic interactions. If a `kspace_style ` is declared, an error will be issued. Related commands diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 80a121c035..d7cbf0362a 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -1560,7 +1560,7 @@ void FixWallGran::unpack_restart(int nlocal, int nth) // skip to Nth set of extra values // unpack the Nth first values this way because other fixes pack them - + int m = 0; for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); m++; diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index 6953165af6..543df02fc2 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -498,7 +498,7 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth) // skip to Nth set of extra values // unpack the Nth first values this way because other fixes pack them - + int m = 0; for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); m++; diff --git a/src/KOKKOS/gridcomm_kokkos.cpp b/src/KOKKOS/gridcomm_kokkos.cpp index 21f2bb915f..c4e4a55da2 100644 --- a/src/KOKKOS/gridcomm_kokkos.cpp +++ b/src/KOKKOS/gridcomm_kokkos.cpp @@ -48,9 +48,9 @@ enum{REGULAR,TILED}; template GridCommKokkos::GridCommKokkos(LAMMPS *lmp, MPI_Comm gcomm, - int gnx, int gny, int gnz, - int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, - int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi) + int gnx, int gny, int gnz, + int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, + int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi) : GridComm(lmp, gcomm, gnx, gny, gnz, ixlo,ixhi, iylo, iyhi, izlo, izhi, @@ -72,9 +72,9 @@ GridCommKokkos::GridCommKokkos(LAMMPS *lmp, MPI_Comm gcomm, template GridCommKokkos::GridCommKokkos(LAMMPS *lmp, MPI_Comm gcomm, int /*flag*/, - int gnx, int gny, int gnz, - int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, - int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi, + int gnx, int gny, int gnz, + int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, + int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi, int /*exlo*/, int /*exhi*/, int /*eylo*/, int /*eyhi*/, int /*ezlo*/, int /*ezhi*/) : GridComm(lmp, gcomm, gnx, gny, gnz, @@ -89,14 +89,14 @@ template GridCommKokkos::~GridCommKokkos() { // regular comm data struct - + for (int i = 0; i < nswap; i++) { swap[i].packlist = NULL; swap[i].unpacklist = NULL; } // tiled comm data structs - + for (int i = 0; i < nsend; i++) send[i].packlist = NULL; @@ -163,7 +163,7 @@ void GridCommKokkos::setup_regular(int &nbuf1, int &nbuf2) // setup swaps = exchange of grid data with one of 6 neighobr procs // can be more than one in a direction if ghost region extends beyond neigh proc // all procs have same swap count, but swapsize npack/nunpack can be empty - + nswap = 0; // send own grid pts to -x processor, recv ghost grid pts from +x processor @@ -420,7 +420,7 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) // access CommTiled to get cut dimension // cut = this proc's inlo in that dim // dim is -1 for proc 0, but never accessed - + rcbinfo = (RCBinfo *) memory->smalloc(nprocs*sizeof(RCBinfo),"GridComm:rcbinfo"); RCBinfo rcbone; @@ -435,14 +435,14 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) // accounts for crossings of periodic boundaries // noverlap = # of overlaps, including self // overlap = vector of overlap info using Overlap data struct - + ghostbox[0] = outxlo; ghostbox[1] = outxhi; ghostbox[2] = outylo; ghostbox[3] = outyhi; ghostbox[4] = outzlo; ghostbox[5] = outzhi; - + pbc[0] = pbc[1] = pbc[2] = 0; memory->create(overlap_procs,nprocs,"GridComm:overlap_procs"); @@ -459,10 +459,10 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) memory->create(proclist,noverlap,"GridComm:proclist"); srequest = (Request *) memory->smalloc(noverlap*sizeof(Request),"GridComm:srequest"); - + int nsend_request = 0; ncopy = 0; - + for (m = 0; m < noverlap; m++) { if (overlap[m].proc == me) ncopy++; else { @@ -470,7 +470,7 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) srequest[nsend_request].sender = me; srequest[nsend_request].index = m; for (i = 0; i < 6; i++) - srequest[nsend_request].box[i] = overlap[m].box[i]; + srequest[nsend_request].box[i] = overlap[m].box[i]; nsend_request++; } } @@ -481,7 +481,7 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) (Request *) memory->smalloc(nrecv_request*sizeof(Request),"GridComm:rrequest"); irregular->exchange_data((char *) srequest,sizeof(Request),(char *) rrequest); irregular->destroy_data(); - + // compute overlaps between received ghost boxes and my owned box // overlap box used to setup my Send data struct and respond to requests @@ -515,7 +515,7 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) } nsend = nrecv_request; - + // reply to each Request message with a Response message // content: index for the overlap on requestor, overlap box on my owned grid @@ -530,13 +530,13 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) // process received responses // box used to setup my Recv data struct after unwrapping via PBC // adjacent = 0 if any box of ghost cells does not adjoin my owned cells - + recv = (Recv *) memory->smalloc(nrecv_response*sizeof(Recv),"GridComm:recv"); k_recv_unpacklist = DAT::tdual_int_2d("GridComm:recv_unpacklist",nrecv_response,k_recv_unpacklist.extent(1)); adjacent = 1; - + for (i = 0; i < nrecv_response; i++) { m = rresponse[i].index; recv[i].proc = overlap[m].proc; @@ -547,21 +547,21 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) zlo = rresponse[i].box[4] + overlap[m].pbc[2] * nz; zhi = rresponse[i].box[5] + overlap[m].pbc[2] * nz; recv[i].nunpack = indices(k_recv_unpacklist,i,xlo,xhi,ylo,yhi,zlo,zhi); - + if (xlo != inxhi+1 && xhi != inxlo-1 && - ylo != inyhi+1 && yhi != inylo-1 && - zlo != inzhi+1 && zhi != inzlo-1) adjacent = 0; + ylo != inyhi+1 && yhi != inylo-1 && + zlo != inzhi+1 && zhi != inzlo-1) adjacent = 0; } nrecv = nrecv_response; // create Copy data struct from overlaps with self - + copy = (Copy *) memory->smalloc(ncopy*sizeof(Copy),"GridComm:copy"); k_copy_packlist = DAT::tdual_int_2d("GridComm:copy_packlist",ncopy,k_copy_packlist.extent(1)); k_copy_unpacklist = DAT::tdual_int_2d("GridComm:copy_unpacklist",ncopy,k_copy_unpacklist.extent(1)); - + ncopy = 0; for (m = 0; m < noverlap; m++) { if (overlap[m].proc != me) continue; @@ -600,7 +600,7 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) int nrequest = MAX(nsend,nrecv); requests = new MPI_Request[nrequest]; - + // clean-up memory->sfree(rcbinfo); @@ -614,7 +614,7 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) // nbuf1 = largest pack or unpack in any Send or Recv or Copy // nbuf2 = larget of sum of all packs or unpacks in Send or Recv - + nbuf1 = 0; for (m = 0; m < ncopy; m++) { @@ -643,7 +643,7 @@ void GridCommKokkos::setup_tiled(int &nbuf1, int &nbuf2) template void GridCommKokkos::forward_comm_kspace(KSpace *kspace, int nper, int which, - FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) + FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) { if (layout == REGULAR) forward_comm_kspace_regular(kspace,nper,which,k_buf1,k_buf2,datatype); @@ -656,7 +656,7 @@ void GridCommKokkos::forward_comm_kspace(KSpace *kspace, int nper, i template void GridCommKokkos:: forward_comm_kspace_regular(KSpace *kspace, int nper, int which, - FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) + FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) { int m; MPI_Request request; @@ -687,9 +687,9 @@ forward_comm_kspace_regular(KSpace *kspace, int nper, int which, } if (swap[m].nunpack) MPI_Irecv(buf2,nper*swap[m].nunpack,datatype, - swap[m].recvproc,0,gridcomm,&request); + swap[m].recvproc,0,gridcomm,&request); if (swap[m].npack) MPI_Send(buf1,nper*swap[m].npack,datatype, - swap[m].sendproc,0,gridcomm); + swap[m].sendproc,0,gridcomm); if (swap[m].nunpack) MPI_Wait(&request,MPI_STATUS_IGNORE); if (!lmp->kokkos->cuda_aware_flag) { @@ -708,7 +708,7 @@ forward_comm_kspace_regular(KSpace *kspace, int nper, int which, template void GridCommKokkos:: forward_comm_kspace_tiled(KSpace *kspace, int nper, int which, - FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) + FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) { int i,m,offset; @@ -724,11 +724,11 @@ forward_comm_kspace_tiled(KSpace *kspace, int nper, int which, } // post all receives - + for (m = 0; m < nrecv; m++) { offset = nper * recv[m].offset; MPI_Irecv(&buf2[offset],nper*recv[m].nunpack,datatype, - recv[m].proc,0,gridcomm,&requests[m]); + recv[m].proc,0,gridcomm,&requests[m]); } // perform all sends to other procs @@ -753,7 +753,7 @@ forward_comm_kspace_tiled(KSpace *kspace, int nper, int which, } // unpack all received data - + for (i = 0; i < nrecv; i++) { MPI_Waitany(nrecv,requests,&m,MPI_STATUS_IGNORE); @@ -764,7 +764,7 @@ forward_comm_kspace_tiled(KSpace *kspace, int nper, int which, offset = nper * recv[m].offset; kspaceKKBase->unpack_forward_grid_kokkos(which,k_buf2,offset, - recv[m].nunpack,k_recv_unpacklist,m); + recv[m].nunpack,k_recv_unpacklist,m); DeviceType().fence(); } } @@ -776,7 +776,7 @@ forward_comm_kspace_tiled(KSpace *kspace, int nper, int which, template void GridCommKokkos::reverse_comm_kspace(KSpace *kspace, int nper, int which, - FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) + FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) { if (layout == REGULAR) reverse_comm_kspace_regular(kspace,nper,which,k_buf1,k_buf2,datatype); @@ -789,7 +789,7 @@ void GridCommKokkos::reverse_comm_kspace(KSpace *kspace, int nper, i template void GridCommKokkos:: reverse_comm_kspace_regular(KSpace *kspace, int nper, int which, - FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) + FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) { int m; MPI_Request request; @@ -820,9 +820,9 @@ reverse_comm_kspace_regular(KSpace *kspace, int nper, int which, } if (swap[m].npack) MPI_Irecv(buf2,nper*swap[m].npack,datatype, - swap[m].sendproc,0,gridcomm,&request); + swap[m].sendproc,0,gridcomm,&request); if (swap[m].nunpack) MPI_Send(buf1,nper*swap[m].nunpack,datatype, - swap[m].recvproc,0,gridcomm); + swap[m].recvproc,0,gridcomm); if (swap[m].npack) MPI_Wait(&request,MPI_STATUS_IGNORE); @@ -842,7 +842,7 @@ reverse_comm_kspace_regular(KSpace *kspace, int nper, int which, template void GridCommKokkos:: reverse_comm_kspace_tiled(KSpace *kspace, int nper, int which, - FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) + FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) { int i,m,offset; @@ -859,11 +859,11 @@ reverse_comm_kspace_tiled(KSpace *kspace, int nper, int which, } // post all receives - + for (m = 0; m < nsend; m++) { offset = nper * send[m].offset; MPI_Irecv(&buf2[offset],nper*send[m].npack,datatype, - send[m].proc,0,gridcomm,&requests[m]); + send[m].proc,0,gridcomm,&requests[m]); } // perform all sends to other procs @@ -888,7 +888,7 @@ reverse_comm_kspace_tiled(KSpace *kspace, int nper, int which, } // unpack all received data - + for (i = 0; i < nsend; i++) { MPI_Waitany(nsend,requests,&m,MPI_STATUS_IGNORE); @@ -899,7 +899,7 @@ reverse_comm_kspace_tiled(KSpace *kspace, int nper, int which, offset = nper * send[m].offset; kspaceKKBase->unpack_reverse_grid_kokkos(which,k_buf2,offset, - send[m].npack,k_send_packlist,m); + send[m].npack,k_send_packlist,m); DeviceType().fence(); } } @@ -942,7 +942,7 @@ int GridCommKokkos::indices(DAT::tdual_int_2d &k_list, int index, if (k_list.extent(1) < nmax) k_list.resize(k_list.extent(0),nmax); - if (nmax == 0) return 0; + if (nmax == 0) return 0; int nx = (fullxhi-fullxlo+1); int ny = (fullyhi-fullylo+1); diff --git a/src/KOKKOS/gridcomm_kokkos.h b/src/KOKKOS/gridcomm_kokkos.h index 1c82ab18e4..1f93c111ca 100644 --- a/src/KOKKOS/gridcomm_kokkos.h +++ b/src/KOKKOS/gridcomm_kokkos.h @@ -27,17 +27,17 @@ class GridCommKokkos : public GridComm { typedef ArrayTypes AT; typedef FFTArrayTypes FFT_AT; GridCommKokkos(class LAMMPS *, MPI_Comm, int, int, int, - int, int, int, int, int, int, - int, int, int, int, int, int); + int, int, int, int, int, int, + int, int, int, int, int, int); GridCommKokkos(class LAMMPS *, MPI_Comm, int, int, int, int, - int, int, int, int, int, int, - int, int, int, int, int, int, - int, int, int, int, int, int); + int, int, int, int, int, int, + int, int, int, int, int, int, + int, int, int, int, int, int); virtual ~GridCommKokkos(); void forward_comm_kspace(class KSpace *, int, int, - FFT_DAT::tdual_FFT_SCALAR_1d &, FFT_DAT::tdual_FFT_SCALAR_1d &, MPI_Datatype); + FFT_DAT::tdual_FFT_SCALAR_1d &, FFT_DAT::tdual_FFT_SCALAR_1d &, MPI_Datatype); void reverse_comm_kspace(class KSpace *, int, int, - FFT_DAT::tdual_FFT_SCALAR_1d &, FFT_DAT::tdual_FFT_SCALAR_1d &, MPI_Datatype); + FFT_DAT::tdual_FFT_SCALAR_1d &, FFT_DAT::tdual_FFT_SCALAR_1d &, MPI_Datatype); private: DAT::tdual_int_2d k_swap_packlist; @@ -67,7 +67,7 @@ class GridCommKokkos : public GridComm { FFT_DAT::tdual_FFT_SCALAR_1d &, FFT_DAT::tdual_FFT_SCALAR_1d &, MPI_Datatype); void grow_swap(); - + int indices(DAT::tdual_int_2d &, int, int, int, int, int, int, int); }; diff --git a/src/KOKKOS/pair_snap_kokkos_impl.h b/src/KOKKOS/pair_snap_kokkos_impl.h index e2a66ee28c..fbdb23a79b 100644 --- a/src/KOKKOS/pair_snap_kokkos_impl.h +++ b/src/KOKKOS/pair_snap_kokkos_impl.h @@ -718,7 +718,7 @@ void PairSNAPKokkos::operator() (TagPairSNAPTransformUi,const int ia const int iatom = iatom_mod + iatom_div * 32; if (iatom >= chunk_size) return; - if (j > twojmax) return; + if (j > twojmax) return; int elem_count = chemflag ? nelements : 1; @@ -739,7 +739,7 @@ void PairSNAPKokkos::operator() (TagPairSNAPTransformUi,const int ia // Store my_sna.ulisttot_pack(iatom_mod, idxu, ielem, iatom_div) = { utot_re, utot_im }; - + // Also zero yi my_sna.ylist_pack_re(iatom_mod, idxu_half, ielem, iatom_div) = 0.; my_sna.ylist_pack_im(iatom_mod, idxu_half, ielem, iatom_div) = 0.; @@ -944,7 +944,7 @@ void PairSNAPKokkos::operator() (TagPairSNAPTransformUiCPU, const in if (iatom >= chunk_size) return; - if (j > twojmax) return; + if (j > twojmax) return; int elem_count = chemflag ? nelements : 1; diff --git a/src/KOKKOS/sna_kokkos.h b/src/KOKKOS/sna_kokkos.h index 765ec1a05a..2c712a1c84 100644 --- a/src/KOKKOS/sna_kokkos.h +++ b/src/KOKKOS/sna_kokkos.h @@ -55,7 +55,7 @@ public: typedef Kokkos::View t_sna_3c3; typedef Kokkos::View t_sna_5c; - typedef Kokkos::View t_sna_2ckp; + typedef Kokkos::View t_sna_2ckp; inline SNAKokkos() {}; diff --git a/src/KOKKOS/sna_kokkos_impl.h b/src/KOKKOS/sna_kokkos_impl.h index 817617090d..0a7dae04a4 100644 --- a/src/KOKKOS/sna_kokkos_impl.h +++ b/src/KOKKOS/sna_kokkos_impl.h @@ -120,7 +120,7 @@ void SNAKokkos::build_indexlist() idxu_max = idxu_count; Kokkos::deep_copy(idxu_block,h_idxu_block); - // index list for half uarray + // index list for half uarray idxu_half_block = Kokkos::View("SNAKokkos::idxu_half_block",jdim); auto h_idxu_half_block = Kokkos::create_mirror_view(idxu_half_block); @@ -316,7 +316,7 @@ void SNAKokkos::grow_rij(int newnatom, int newnmax) /* ---------------------------------------------------------------------- Precompute the Cayley-Klein parameters and the derivatives thereof. - This routine better exploits parallelism than the GPU ComputeUi and + This routine better exploits parallelism than the GPU ComputeUi and ComputeFusedDeidrj, which are one warp per atom-neighbor pair. ------------------------------------------------------------------------- */ @@ -339,7 +339,7 @@ void SNAKokkos::compute_cayley_klein(const int& iatom, const int& jn const double dz0dr = z0 / r - (r*rscale0) * (rsq + z0 * z0) / rsq; //const double wj_local = wj(iatom, jnbor); - double sfac, dsfac; + double sfac, dsfac; compute_s_dsfac(r, rcut, sfac, dsfac); sfac *= wj_local; dsfac *= wj_local; @@ -639,7 +639,7 @@ void SNAKokkos::compute_zi(const int& iatom_mod, const int& jjz, con /* ---------------------------------------------------------------------- compute Bi by summing conj(Ui)*Zi AoSoA data layout to take advantage of coalescing, avoiding warp - divergence. + divergence. ------------------------------------------------------------------------- */ template @@ -1406,7 +1406,7 @@ void SNAKokkos::compute_deidrj_cpu(const typename Kokkos::TeamPolicy for(int ma = 0; ma <= j; ma++) { sum_tmp.x += dulist(jju_cache,iatom,jnbor,0).re * ylist(jju_half,jelem,iatom).re + dulist(jju_cache,iatom,jnbor,0).im * ylist(jju_half,jelem,iatom).im; - sum_tmp.y += dulist(jju_cache,iatom,jnbor,1).re * ylist(jju_half,jelem,iatom).re + + sum_tmp.y += dulist(jju_cache,iatom,jnbor,1).re * ylist(jju_half,jelem,iatom).re + dulist(jju_cache,iatom,jnbor,1).im * ylist(jju_half,jelem,iatom).im; sum_tmp.z += dulist(jju_cache,iatom,jnbor,2).re * ylist(jju_half,jelem,iatom).re + dulist(jju_cache,iatom,jnbor,2).im * ylist(jju_half,jelem,iatom).im; @@ -1421,9 +1421,9 @@ void SNAKokkos::compute_deidrj_cpu(const typename Kokkos::TeamPolicy for(int ma = 0; ma < mb; ma++) { sum_tmp.x += dulist(jju_cache,iatom,jnbor,0).re * ylist(jju_half,jelem,iatom).re + dulist(jju_cache,iatom,jnbor,0).im * ylist(jju_half,jelem,iatom).im; - sum_tmp.y += dulist(jju_cache,iatom,jnbor,1).re * ylist(jju_half,jelem,iatom).re + + sum_tmp.y += dulist(jju_cache,iatom,jnbor,1).re * ylist(jju_half,jelem,iatom).re + dulist(jju_cache,iatom,jnbor,1).im * ylist(jju_half,jelem,iatom).im; - sum_tmp.z += dulist(jju_cache,iatom,jnbor,2).re * ylist(jju_half,jelem,iatom).re + + sum_tmp.z += dulist(jju_cache,iatom,jnbor,2).re * ylist(jju_half,jelem,iatom).re + dulist(jju_cache,iatom,jnbor,2).im * ylist(jju_half,jelem,iatom).im; jju_half++; jju_cache++; } @@ -1431,9 +1431,9 @@ void SNAKokkos::compute_deidrj_cpu(const typename Kokkos::TeamPolicy //int ma = mb; sum_tmp.x += (dulist(jju_cache,iatom,jnbor,0).re * ylist(jju_half,jelem,iatom).re + dulist(jju_cache,iatom,jnbor,0).im * ylist(jju_half,jelem,iatom).im)*0.5; - sum_tmp.y += (dulist(jju_cache,iatom,jnbor,1).re * ylist(jju_half,jelem,iatom).re + + sum_tmp.y += (dulist(jju_cache,iatom,jnbor,1).re * ylist(jju_half,jelem,iatom).re + dulist(jju_cache,iatom,jnbor,1).im * ylist(jju_half,jelem,iatom).im)*0.5; - sum_tmp.z += (dulist(jju_cache,iatom,jnbor,2).re * ylist(jju_half,jelem,iatom).re + + sum_tmp.z += (dulist(jju_cache,iatom,jnbor,2).re * ylist(jju_half,jelem,iatom).re + dulist(jju_cache,iatom,jnbor,2).im * ylist(jju_half,jelem,iatom).im)*0.5; } // end if jeven @@ -2162,7 +2162,7 @@ double SNAKokkos::memory_usage() #ifdef LMP_KOKKOS_GPU if (std::is_same::value) { - + auto natom_pad = (natom+32-1)/32; bytes += natom * idxu_half_max * nelements * sizeof(double); // ulisttot_re diff --git a/src/KSPACE/fft3d.cpp b/src/KSPACE/fft3d.cpp index 0b4bb85032..7c555e99b5 100644 --- a/src/KSPACE/fft3d.cpp +++ b/src/KSPACE/fft3d.cpp @@ -198,7 +198,7 @@ void fft_3d(FFT_DATA *in, FFT_DATA *out, int flag, struct fft_plan_3d *plan) (FFT_SCALAR *) plan->scratch, plan->post_plan); // scaling if required - + if (flag == -1 && plan->scaled) { norm = plan->norm; num = plan->normnum; diff --git a/src/KSPACE/gridcomm.cpp b/src/KSPACE/gridcomm.cpp index 27952879ce..4024670b9f 100644 --- a/src/KSPACE/gridcomm.cpp +++ b/src/KSPACE/gridcomm.cpp @@ -46,9 +46,9 @@ enum{REGULAR,TILED}; ------------------------------------------------------------------------- */ GridComm::GridComm(LAMMPS *lmp, MPI_Comm gcomm, - int gnx, int gny, int gnz, - int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, - int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi) + int gnx, int gny, int gnz, + int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, + int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi) : Pointers(lmp) { if (comm->layout == Comm::LAYOUT_TILED) layout = TILED; @@ -57,18 +57,18 @@ GridComm::GridComm(LAMMPS *lmp, MPI_Comm gcomm, if (layout == REGULAR) { int (*procneigh)[2] = comm->procneigh; initialize(gcomm,gnx,gny,gnz, - ixlo,ixhi,iylo,iyhi,izlo,izhi, - oxlo,oxhi,oylo,oyhi,ozlo,ozhi, - oxlo,oxhi,oylo,oyhi,ozlo,ozhi, - procneigh[0][0],procneigh[0][1], - procneigh[1][0],procneigh[1][1], - procneigh[2][0],procneigh[2][1]); + ixlo,ixhi,iylo,iyhi,izlo,izhi, + oxlo,oxhi,oylo,oyhi,ozlo,ozhi, + oxlo,oxhi,oylo,oyhi,ozlo,ozhi, + procneigh[0][0],procneigh[0][1], + procneigh[1][0],procneigh[1][1], + procneigh[2][0],procneigh[2][1]); } else { initialize(gcomm,gnx,gny,gnz, - ixlo,ixhi,iylo,iyhi,izlo,izhi, - oxlo,oxhi,oylo,oyhi,ozlo,ozhi, - oxlo,oxhi,oylo,oyhi,ozlo,ozhi, - 0,0,0,0,0,0); + ixlo,ixhi,iylo,iyhi,izlo,izhi, + oxlo,oxhi,oylo,oyhi,ozlo,ozhi, + oxlo,oxhi,oylo,oyhi,ozlo,ozhi, + 0,0,0,0,0,0); } } @@ -85,10 +85,10 @@ GridComm::GridComm(LAMMPS *lmp, MPI_Comm gcomm, ------------------------------------------------------------------------- */ GridComm::GridComm(LAMMPS *lmp, MPI_Comm gcomm, int flag, - int gnx, int gny, int gnz, - int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, - int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi, - int exlo, int exhi, int eylo, int eyhi, int ezlo, int ezhi) + int gnx, int gny, int gnz, + int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, + int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi, + int exlo, int exhi, int eylo, int eyhi, int ezlo, int ezhi) : Pointers(lmp) { if (comm->layout == Comm::LAYOUT_TILED) layout = TILED; @@ -99,27 +99,27 @@ GridComm::GridComm(LAMMPS *lmp, MPI_Comm gcomm, int flag, // this assumes gcomm = world int (*procneigh)[2] = comm->procneigh; initialize(gcomm,gnx,gny,gnz, - ixlo,ixhi,iylo,iyhi,izlo,izhi, - oxlo,oxhi,oylo,oyhi,ozlo,ozhi, - exlo,exhi,eylo,eyhi,ezlo,ezhi, - procneigh[0][0],procneigh[0][1], - procneigh[1][0],procneigh[1][1], - procneigh[2][0],procneigh[2][1]); + ixlo,ixhi,iylo,iyhi,izlo,izhi, + oxlo,oxhi,oylo,oyhi,ozlo,ozhi, + exlo,exhi,eylo,eyhi,ezlo,ezhi, + procneigh[0][0],procneigh[0][1], + procneigh[1][0],procneigh[1][1], + procneigh[2][0],procneigh[2][1]); } else { initialize(gcomm,gnx,gny,gnz, - ixlo,ixhi,iylo,iyhi,izlo,izhi, - oxlo,oxhi,oylo,oyhi,ozlo,ozhi, - exlo,exhi,eylo,eyhi,ezlo,ezhi, - 0,0,0,0,0,0); + ixlo,ixhi,iylo,iyhi,izlo,izhi, + oxlo,oxhi,oylo,oyhi,ozlo,ozhi, + exlo,exhi,eylo,eyhi,ezlo,ezhi, + 0,0,0,0,0,0); } - + } else if (flag == 2) { if (layout == REGULAR) { initialize(gcomm,gnx,gny,gnz, - ixlo,ixhi,iylo,iyhi,izlo,izhi, - oxlo,oxhi,oylo,oyhi,ozlo,ozhi, - oxlo,oxhi,oylo,oyhi,ozlo,ozhi, - exlo,exhi,eylo,eyhi,ezlo,ezhi); + ixlo,ixhi,iylo,iyhi,izlo,izhi, + oxlo,oxhi,oylo,oyhi,ozlo,ozhi, + oxlo,oxhi,oylo,oyhi,ozlo,ozhi, + exlo,exhi,eylo,eyhi,ezlo,ezhi); } else { error->all(FLERR,"GridComm does not support tiled layout with neighbor procs"); } @@ -131,7 +131,7 @@ GridComm::GridComm(LAMMPS *lmp, MPI_Comm gcomm, int flag, GridComm::~GridComm() { // regular comm data struct - + for (int i = 0; i < nswap; i++) { memory->destroy(swap[i].packlist); memory->destroy(swap[i].unpacklist); @@ -139,7 +139,7 @@ GridComm::~GridComm() memory->sfree(swap); // tiled comm data structs - + for (int i = 0; i < nsend; i++) memory->destroy(send[i].packlist); memory->sfree(send); @@ -162,11 +162,11 @@ GridComm::~GridComm() ------------------------------------------------------------------------- */ void GridComm::initialize(MPI_Comm gcomm, - int gnx, int gny, int gnz, - int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, - int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi, - int fxlo, int fxhi, int fylo, int fyhi, int fzlo, int fzhi, - int pxlo, int pxhi, int pylo, int pyhi, int pzlo, int pzhi) + int gnx, int gny, int gnz, + int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, + int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi, + int fxlo, int fxhi, int fylo, int fyhi, int fzlo, int fzhi, + int pxlo, int pxhi, int pylo, int pyhi, int pzlo, int pzhi) { gridcomm = gcomm; MPI_Comm_rank(gridcomm,&me); @@ -175,7 +175,7 @@ void GridComm::initialize(MPI_Comm gcomm, nx = gnx; ny = gny; nz = gnz; - + inxlo = ixlo; inxhi = ixhi; inylo = iylo; @@ -209,7 +209,7 @@ void GridComm::initialize(MPI_Comm gcomm, } // internal data initializations - + nswap = maxswap = 0; swap = NULL; @@ -287,7 +287,7 @@ void GridComm::setup_regular(int &nbuf1, int &nbuf2) // setup swaps = exchange of grid data with one of 6 neighobr procs // can be more than one in a direction if ghost region extends beyond neigh proc // all procs have same swap count, but swapsize npack/nunpack can be empty - + nswap = 0; // send own grid pts to -x processor, recv ghost grid pts from +x processor @@ -548,7 +548,7 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) // access CommTiled to get cut dimension // cut = this proc's inlo in that dim // dim is -1 for proc 0, but never accessed - + rcbinfo = (RCBinfo *) memory->smalloc(nprocs*sizeof(RCBinfo),"GridComm:rcbinfo"); RCBinfo rcbone; @@ -563,14 +563,14 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) // accounts for crossings of periodic boundaries // noverlap = # of overlaps, including self // overlap = vector of overlap info using Overlap data struct - + ghostbox[0] = outxlo; ghostbox[1] = outxhi; ghostbox[2] = outylo; ghostbox[3] = outyhi; ghostbox[4] = outzlo; ghostbox[5] = outzhi; - + pbc[0] = pbc[1] = pbc[2] = 0; memory->create(overlap_procs,nprocs,"GridComm:overlap_procs"); @@ -587,10 +587,10 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) memory->create(proclist,noverlap,"GridComm:proclist"); srequest = (Request *) memory->smalloc(noverlap*sizeof(Request),"GridComm:srequest"); - + int nsend_request = 0; ncopy = 0; - + for (m = 0; m < noverlap; m++) { if (overlap[m].proc == me) ncopy++; else { @@ -598,7 +598,7 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) srequest[nsend_request].sender = me; srequest[nsend_request].index = m; for (i = 0; i < 6; i++) - srequest[nsend_request].box[i] = overlap[m].box[i]; + srequest[nsend_request].box[i] = overlap[m].box[i]; nsend_request++; } } @@ -609,7 +609,7 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) (Request *) memory->smalloc(nrecv_request*sizeof(Request),"GridComm:rrequest"); irregular->exchange_data((char *) srequest,sizeof(Request),(char *) rrequest); irregular->destroy_data(); - + // compute overlaps between received ghost boxes and my owned box // overlap box used to setup my Send data struct and respond to requests @@ -640,7 +640,7 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) } nsend = nrecv_request; - + // reply to each Request message with a Response message // content: index for the overlap on requestor, overlap box on my owned grid @@ -655,10 +655,10 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) // process received responses // box used to setup my Recv data struct after unwrapping via PBC // adjacent = 0 if any box of ghost cells does not adjoin my owned cells - + recv = (Recv *) memory->smalloc(nrecv_response*sizeof(Recv),"GridComm:recv"); adjacent = 1; - + for (i = 0; i < nrecv_response; i++) { m = rresponse[i].index; recv[i].proc = overlap[m].proc; @@ -669,18 +669,18 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) zlo = rresponse[i].box[4] + overlap[m].pbc[2] * nz; zhi = rresponse[i].box[5] + overlap[m].pbc[2] * nz; recv[i].nunpack = indices(recv[i].unpacklist,xlo,xhi,ylo,yhi,zlo,zhi); - + if (xlo != inxhi+1 && xhi != inxlo-1 && - ylo != inyhi+1 && yhi != inylo-1 && - zlo != inzhi+1 && zhi != inzlo-1) adjacent = 0; + ylo != inyhi+1 && yhi != inylo-1 && + zlo != inzhi+1 && zhi != inzlo-1) adjacent = 0; } nrecv = nrecv_response; // create Copy data struct from overlaps with self - + copy = (Copy *) memory->smalloc(ncopy*sizeof(Copy),"GridComm:copy"); - + ncopy = 0; for (m = 0; m < noverlap; m++) { if (overlap[m].proc != me) continue; @@ -719,7 +719,7 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) int nrequest = MAX(nsend,nrecv); requests = new MPI_Request[nrequest]; - + // clean-up memory->sfree(rcbinfo); @@ -733,7 +733,7 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) // nbuf1 = largest pack or unpack in any Send or Recv or Copy // nbuf2 = larget of sum of all packs or unpacks in Send or Recv - + nbuf1 = 0; for (m = 0; m < ncopy; m++) { @@ -769,9 +769,9 @@ void GridComm::setup_tiled(int &nbuf1, int &nbuf2) void GridComm::ghost_box_drop(int *box, int *pbc) { int i,m; - + // newbox12 and newpbc are initially copies of caller box and pbc - + int newbox1[6],newbox2[6],newpbc[3]; for (i = 0; i < 6; i++) newbox1[i] = newbox2[i] = box[i]; @@ -780,9 +780,9 @@ void GridComm::ghost_box_drop(int *box, int *pbc) // 6 if tests to see if box needs to be split across a periodic boundary // newbox1 and 2 = new split boxes, newpbc increments current pbc // final else is no split - + int splitflag = 1; - + if (box[0] < 0) { newbox1[0] = 0; newbox2[0] = box[0] + nx; @@ -801,7 +801,7 @@ void GridComm::ghost_box_drop(int *box, int *pbc) } else if (box[3] >= ny) { newbox1[3] = ny - 1; newbox2[2] = 0; - newbox2[3] = box[3] - ny; + newbox2[3] = box[3] - ny; newpbc[1]++; } else if (box[4] < 0) { newbox1[4] = 0; @@ -819,7 +819,7 @@ void GridComm::ghost_box_drop(int *box, int *pbc) // returns proc_overlap = list of proc IDs it overlaps // skip self overlap if no crossing of periodic boundaries // do not skip self if overlap is in another periodic image - + } else { splitflag = 0; int np = 0; @@ -827,7 +827,7 @@ void GridComm::ghost_box_drop(int *box, int *pbc) for (m = 0; m < np; m++) { if (noverlap == maxoverlap) grow_overlap(); if (overlap_procs[m] == me && - pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 0) continue; + pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 0) continue; overlap[noverlap].proc = overlap_procs[m]; for (i = 0; i < 6; i++) overlap[noverlap].box[i] = box[i]; for (i = 0; i < 3; i++) overlap[noverlap].pbc[i] = pbc[i]; @@ -836,7 +836,7 @@ void GridComm::ghost_box_drop(int *box, int *pbc) } // recurse with 2 split boxes - + if (splitflag) { ghost_box_drop(newbox1,pbc); ghost_box_drop(newbox2,newpbc); @@ -852,7 +852,7 @@ void GridComm::ghost_box_drop(int *box, int *pbc) ------------------------------------------------------------------------- */ void GridComm::box_drop_grid(int *box, int proclower, int procupper, - int &np, int *plist) + int &np, int *plist) { // end recursion when partition is a single proc // add proclower to plist @@ -926,7 +926,7 @@ int GridComm::ghost_adjacent_tiled() ------------------------------------------------------------------------- */ void GridComm::forward_comm_kspace(KSpace *kspace, int nper, int nbyte, int which, - void *buf1, void *buf2, MPI_Datatype datatype) + void *buf1, void *buf2, MPI_Datatype datatype) { if (layout == REGULAR) forward_comm_kspace_regular(kspace,nper,nbyte,which,buf1,buf2,datatype); @@ -940,7 +940,7 @@ void GridComm::forward_comm_kspace(KSpace *kspace, int nper, int nbyte, int whic void GridComm:: forward_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which, - void *buf1, void *buf2, MPI_Datatype datatype) + void *buf1, void *buf2, MPI_Datatype datatype) { int m; MPI_Request request; @@ -953,9 +953,9 @@ forward_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which, if (swap[m].sendproc != me) { if (swap[m].nunpack) MPI_Irecv(buf2,nper*swap[m].nunpack,datatype, - swap[m].recvproc,0,gridcomm,&request); + swap[m].recvproc,0,gridcomm,&request); if (swap[m].npack) MPI_Send(buf1,nper*swap[m].npack,datatype, - swap[m].sendproc,0,gridcomm); + swap[m].sendproc,0,gridcomm); if (swap[m].nunpack) MPI_Wait(&request,MPI_STATUS_IGNORE); } @@ -969,18 +969,18 @@ forward_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which, void GridComm:: forward_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which, - void *buf1, void *vbuf2, MPI_Datatype datatype) + void *buf1, void *vbuf2, MPI_Datatype datatype) { int i,m,offset; char *buf2 = (char *) vbuf2; - + // post all receives - + for (m = 0; m < nrecv; m++) { offset = nper * recv[m].offset * nbyte; MPI_Irecv((void *) &buf2[offset],nper*recv[m].nunpack,datatype, - recv[m].proc,0,gridcomm,&requests[m]); + recv[m].proc,0,gridcomm,&requests[m]); } // perform all sends to other procs @@ -998,12 +998,12 @@ forward_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which, } // unpack all received data - + for (i = 0; i < nrecv; i++) { MPI_Waitany(nrecv,requests,&m,MPI_STATUS_IGNORE); offset = nper * recv[m].offset * nbyte; kspace->unpack_forward_grid(which,(void *) &buf2[offset], - recv[m].nunpack,recv[m].unpacklist); + recv[m].nunpack,recv[m].unpacklist); } } @@ -1012,7 +1012,7 @@ forward_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which, ------------------------------------------------------------------------- */ void GridComm::reverse_comm_kspace(KSpace *kspace, int nper, int nbyte, int which, - void *buf1, void *buf2, MPI_Datatype datatype) + void *buf1, void *buf2, MPI_Datatype datatype) { if (layout == REGULAR) reverse_comm_kspace_regular(kspace,nper,nbyte,which,buf1,buf2,datatype); @@ -1026,7 +1026,7 @@ void GridComm::reverse_comm_kspace(KSpace *kspace, int nper, int nbyte, int whic void GridComm:: reverse_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which, - void *buf1, void *buf2, MPI_Datatype datatype) + void *buf1, void *buf2, MPI_Datatype datatype) { int m; MPI_Request request; @@ -1039,9 +1039,9 @@ reverse_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which, if (swap[m].recvproc != me) { if (swap[m].npack) MPI_Irecv(buf2,nper*swap[m].npack,datatype, - swap[m].sendproc,0,gridcomm,&request); + swap[m].sendproc,0,gridcomm,&request); if (swap[m].nunpack) MPI_Send(buf1,nper*swap[m].nunpack,datatype, - swap[m].recvproc,0,gridcomm); + swap[m].recvproc,0,gridcomm); if (swap[m].npack) MPI_Wait(&request,MPI_STATUS_IGNORE); } @@ -1055,18 +1055,18 @@ reverse_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which, void GridComm:: reverse_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which, - void *buf1, void *vbuf2, MPI_Datatype datatype) + void *buf1, void *vbuf2, MPI_Datatype datatype) { int i,m,offset; char *buf2 = (char *) vbuf2; // post all receives - + for (m = 0; m < nsend; m++) { offset = nper * send[m].offset * nbyte; MPI_Irecv((void *) &buf2[offset],nper*send[m].npack,datatype, - send[m].proc,0,gridcomm,&requests[m]); + send[m].proc,0,gridcomm,&requests[m]); } // perform all sends to other procs @@ -1084,12 +1084,12 @@ reverse_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which, } // unpack all received data - + for (i = 0; i < nsend; i++) { MPI_Waitany(nsend,requests,&m,MPI_STATUS_IGNORE); offset = nper * send[m].offset * nbyte; kspace->unpack_reverse_grid(which,(void *) &buf2[offset], - send[m].npack,send[m].packlist); + send[m].npack,send[m].packlist); } } @@ -1125,7 +1125,7 @@ void GridComm::grow_overlap() /* ---------------------------------------------------------------------- create 1d list of offsets into 3d array section (xlo:xhi,ylo:yhi,zlo:zhi) - assume 3d array is allocated as + assume 3d array is allocated as (fullxlo:fullxhi,fullylo:fullyhi,fullzlo:fullzhi) ------------------------------------------------------------------------- */ diff --git a/src/KSPACE/gridcomm.h b/src/KSPACE/gridcomm.h index 941b84651f..97c914999f 100644 --- a/src/KSPACE/gridcomm.h +++ b/src/KSPACE/gridcomm.h @@ -21,12 +21,12 @@ namespace LAMMPS_NS { class GridComm : protected Pointers { public: GridComm(class LAMMPS *, MPI_Comm, int, int, int, - int, int, int, int, int, int, - int, int, int, int, int, int); + int, int, int, int, int, int, + int, int, int, int, int, int); GridComm(class LAMMPS *, MPI_Comm, int, int, int, int, - int, int, int, int, int, int, - int, int, int, int, int, int, - int, int, int, int, int, int); + int, int, int, int, int, int, + int, int, int, int, int, int, + int, int, int, int, int, int); virtual ~GridComm(); void setup(int &, int &); int ghost_adjacent(); @@ -46,7 +46,7 @@ class GridComm : protected Pointers { int nx,ny,nz; // size of global grid in all 3 dims int inxlo,inxhi; // inclusive extent of my grid chunk int inylo,inyhi; // 0 <= in <= N-1 - int inzlo,inzhi; + int inzlo,inzhi; int outxlo,outxhi; // inclusive extent of my grid chunk plus int outylo,outyhi; // ghost cells in all 6 directions int outzlo,outzhi; // lo indices can be < 0, hi indices can be >= N @@ -61,13 +61,13 @@ class GridComm : protected Pointers { int procxlo,procxhi; // 6 neighbor procs that adjoin me int procylo,procyhi; // not used for comm_style = tiled int proczlo,proczhi; - + int ghostxlo,ghostxhi; // # of my owned grid planes needed int ghostylo,ghostyhi; // by neighobr procs in each dir as their ghost planes int ghostzlo,ghostzhi; // swap = exchange of owned and ghost grid cells between 2 procs, including self - + struct Swap { int sendproc; // proc to send to for forward comm int recvproc; // proc to recv from for forward comm @@ -89,17 +89,17 @@ class GridComm : protected Pointers { // RCB tree of cut info // each proc contributes one value, except proc 0 - + struct RCBinfo { int dim; // 0,1,2 = which dim the cut is in int cut; // grid index of lowest cell in upper half of cut }; RCBinfo *rcbinfo; - + // overlap = a proc whose owned cells overlap with my extended ghost box // includes overlaps across periodic boundaries, can also be self - + struct Overlap { int proc; // proc whose owned cells overlap my ghost cells int box[6]; // box that overlaps otherproc's owned cells @@ -110,9 +110,9 @@ class GridComm : protected Pointers { int noverlap,maxoverlap; Overlap *overlap; - + // request = sent to each proc whose owned cells overlap my ghost cells - + struct Request { int sender; // sending proc int index; // index of overlap on sender @@ -121,9 +121,9 @@ class GridComm : protected Pointers { }; Request *srequest,*rrequest; - + // response = reply from each proc whose owned cells overlap my ghost cells - + struct Response { int index; // index of my overlap for the initial request int box[6]; // box that overlaps responder's owned cells @@ -132,7 +132,7 @@ class GridComm : protected Pointers { }; Response *sresponse,*rresponse; - + // send = proc to send a subset of my owned cells to, for forward comm // for reverse comm, proc I receive ghost overlaps with my owned cells from // offset used in reverse comm to recv a message in middle of a large buffer @@ -147,7 +147,7 @@ class GridComm : protected Pointers { // recv = proc to recv a subset of my ghost cells from, for forward comm // for reverse comm, proc I send a subset of my ghost cells to // offset used in forward comm to recv a message in middle of a large buffer - + struct Recv { int proc; int nunpack; @@ -159,7 +159,7 @@ class GridComm : protected Pointers { // copy = subset of my owned cells to copy into subset of my ghost cells // that describes forward comm, for reverse comm it is the opposite - + struct Copy { int npack; int nunpack; @@ -177,18 +177,18 @@ class GridComm : protected Pointers { // ------------------------------------------- void initialize(MPI_Comm, int, int, int, - int, int, int, int, int, int, - int, int, int, int, int, int, - int, int, int, int, int, int, - int, int, int, int, int, int); + int, int, int, int, int, int, + int, int, int, int, int, int, + int, int, int, int, int, int, + int, int, int, int, int, int); virtual void setup_regular(int &, int &); virtual void setup_tiled(int &, int &); void ghost_box_drop(int *, int *); void box_drop_grid(int *, int, int, int &, int *); - + int ghost_adjacent_regular(); int ghost_adjacent_tiled(); - + void forward_comm_kspace_regular(class KSpace *, int, int, int, void *, void *, MPI_Datatype); void forward_comm_kspace_tiled(class KSpace *, int, int, int, @@ -200,7 +200,7 @@ class GridComm : protected Pointers { virtual void grow_swap(); void grow_overlap(); - + int indices(int *&, int, int, int, int, int, int); }; diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index f2c3f6c820..973302a054 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -467,7 +467,7 @@ void MSM::compute(int eflag, int vflag) current_level = 0; gcall->reverse_comm_kspace(this,1,sizeof(double),REVERSE_RHO, - gcall_buf1,gcall_buf2,MPI_DOUBLE); + gcall_buf1,gcall_buf2,MPI_DOUBLE); // forward communicate charge density values to fill ghost grid points // compute direct sum interaction and then restrict to coarser grid @@ -476,7 +476,7 @@ void MSM::compute(int eflag, int vflag) if (!active_flag[n]) continue; current_level = n; gc[n]->forward_comm_kspace(this,1,sizeof(double),FORWARD_RHO, - gc_buf1[n],gc_buf2[n],MPI_DOUBLE); + gc_buf1[n],gc_buf2[n],MPI_DOUBLE); direct(n); restriction(n); } @@ -488,17 +488,17 @@ void MSM::compute(int eflag, int vflag) if (domain->nonperiodic) { current_level = levels-1; gc[levels-1]-> - forward_comm_kspace(this,1,sizeof(double),FORWARD_RHO, - gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE); + forward_comm_kspace(this,1,sizeof(double),FORWARD_RHO, + gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE); direct_top(levels-1); gc[levels-1]-> - reverse_comm_kspace(this,1,sizeof(double),REVERSE_AD, - gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE); + reverse_comm_kspace(this,1,sizeof(double),REVERSE_AD, + gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE); if (vflag_atom) - gc[levels-1]-> - reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM, - gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE); - + gc[levels-1]-> + reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM, + gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE); + } else { // Here using MPI_Allreduce is cheaper than using commgrid grid_swap_forward(levels-1,qgrid[levels-1]); @@ -506,9 +506,9 @@ void MSM::compute(int eflag, int vflag) grid_swap_reverse(levels-1,egrid[levels-1]); current_level = levels-1; if (vflag_atom) - gc[levels-1]-> - reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM, - gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE); + gc[levels-1]-> + reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM, + gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE); } } @@ -521,13 +521,13 @@ void MSM::compute(int eflag, int vflag) current_level = n; gc[n]->reverse_comm_kspace(this,1,sizeof(double),REVERSE_AD, - gc_buf1[n],gc_buf2[n],MPI_DOUBLE); + gc_buf1[n],gc_buf2[n],MPI_DOUBLE); // extra per-atom virial communication if (vflag_atom) gc[n]->reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM, - gc_buf1[n],gc_buf2[n],MPI_DOUBLE); + gc_buf1[n],gc_buf2[n],MPI_DOUBLE); } // all procs communicate E-field values @@ -535,13 +535,13 @@ void MSM::compute(int eflag, int vflag) current_level = 0; gcall->forward_comm_kspace(this,1,sizeof(double),FORWARD_AD, - gcall_buf1,gcall_buf2,MPI_DOUBLE); + gcall_buf1,gcall_buf2,MPI_DOUBLE); // extra per-atom energy/virial communication if (vflag_atom) gcall->forward_comm_kspace(this,6,sizeof(double),FORWARD_AD_PERATOM, - gcall_buf1,gcall_buf2,MPI_DOUBLE); + gcall_buf1,gcall_buf2,MPI_DOUBLE); // calculate the force on my particles (interpolation) @@ -618,12 +618,12 @@ void MSM::allocate() // commgrid using all processors for finest grid level gcall = new GridComm(lmp,world,1,nx_msm[0],ny_msm[0],nz_msm[0], - nxlo_in[0],nxhi_in[0],nylo_in[0], - nyhi_in[0],nzlo_in[0],nzhi_in[0], - nxlo_out_all,nxhi_out_all,nylo_out_all, - nyhi_out_all,nzlo_out_all,nzhi_out_all, - nxlo_out[0],nxhi_out[0],nylo_out[0], - nyhi_out[0],nzlo_out[0],nzhi_out[0]); + nxlo_in[0],nxhi_in[0],nylo_in[0], + nyhi_in[0],nzlo_in[0],nzhi_in[0], + nxlo_out_all,nxhi_out_all,nylo_out_all, + nyhi_out_all,nzlo_out_all,nzhi_out_all, + nxlo_out[0],nxhi_out[0],nylo_out[0], + nyhi_out[0],nzlo_out[0],nzhi_out[0]); gcall->setup(ngcall_buf1,ngcall_buf2); npergrid = 1; @@ -644,12 +644,12 @@ void MSM::allocate() if (active_flag[n]) { int **procneigh = procneigh_levels[n]; gc[n] = new GridComm(lmp,world_levels[n],2,nx_msm[n],ny_msm[n],nz_msm[n], - nxlo_in[n],nxhi_in[n],nylo_in[n],nyhi_in[n], - nzlo_in[n],nzhi_in[n], - nxlo_out[n],nxhi_out[n],nylo_out[n],nyhi_out[n], - nzlo_out[n],nzhi_out[n], - procneigh[0][0],procneigh[0][1],procneigh[1][0], - procneigh[1][1],procneigh[2][0],procneigh[2][1]); + nxlo_in[n],nxhi_in[n],nylo_in[n],nyhi_in[n], + nzlo_in[n],nzhi_in[n], + nxlo_out[n],nxhi_out[n],nylo_out[n],nyhi_out[n], + nzlo_out[n],nzhi_out[n], + procneigh[0][0],procneigh[0][1],procneigh[1][0], + procneigh[1][1],procneigh[2][0],procneigh[2][1]); gc[n]->setup(ngc_buf1[n],ngc_buf2[n]); npergrid = 1; @@ -677,7 +677,7 @@ void MSM::deallocate() memory->destroy(gcall_buf2); gcall = nullptr; gcall_buf1 = gcall_buf2 = nullptr; - + for (int n=0; ndestroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); @@ -692,10 +692,10 @@ void MSM::deallocate() if (gc) { if (gc[n]) { delete gc[n]; - memory->destroy(gc_buf1[n]); - memory->destroy(gc_buf2[n]); + memory->destroy(gc_buf1[n]); + memory->destroy(gc_buf2[n]); gc[n] = nullptr; - gc_buf1[n] = gc_buf2[n] = nullptr; + gc_buf1[n] = gc_buf2[n] = nullptr; } } } @@ -776,13 +776,13 @@ void MSM::deallocate_peratom() void MSM::allocate_levels() { ngrid = new int[levels]; - + gc = new GridComm*[levels]; gc_buf1 = new double*[levels]; gc_buf2 = new double*[levels]; ngc_buf1 = new int[levels]; ngc_buf2 = new int[levels]; - + memory->create(procneigh_levels,levels,3,2,"msm:procneigh_levels"); world_levels = new MPI_Comm[levels]; active_flag = new int[levels]; @@ -2546,7 +2546,7 @@ void MSM::grid_swap_reverse(int n, double*** &gridn) void MSM::pack_forward_grid(int flag, void *vbuf, int nlist, int *list) { double *buf = (double *) vbuf; - + int n = current_level; int k = 0; @@ -3432,7 +3432,7 @@ double MSM::memory_usage() double bytes = 0; // NOTE: Stan, fill in other memory allocations here - + // all GridComm bufs bytes += (ngcall_buf1 + ngcall_buf2) * npergrid * sizeof(double); diff --git a/src/KSPACE/msm.h b/src/KSPACE/msm.h index 12b0cbf309..a239e6f139 100644 --- a/src/KSPACE/msm.h +++ b/src/KSPACE/msm.h @@ -80,7 +80,7 @@ class MSM : public KSpace { int procgrid[3]; // procs assigned in each dim of 3d grid int myloc[3]; // which proc I am in each dim int ***procneigh_levels; // my 6 neighboring procs, 0/1 = left/right - + class GridComm *gcall; // GridComm class for finest level grid class GridComm **gc; // GridComm classes for each hierarchical level @@ -133,7 +133,7 @@ class MSM : public KSpace { void get_virial_direct_top(int); // grid communication - + void pack_forward_grid(int, void *, int, int *); void unpack_forward_grid(int, void *, int, int *); void pack_reverse_grid(int, void *, int, int *); diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 84b5d9ecb0..e399727001 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -309,8 +309,8 @@ void PPPM::init() if (overlap_allowed) break; gctmp = new GridComm(lmp,world,nx_pppm,ny_pppm,nz_pppm, - nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, - nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); + nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, + nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); int tmp1,tmp2; gctmp->setup(tmp1,tmp2); @@ -641,7 +641,7 @@ void PPPM::compute(int eflag, int vflag) // remap from 3d decomposition to FFT decomposition gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft(); // compute potential gradient on my FFT grid and @@ -656,20 +656,20 @@ void PPPM::compute(int eflag, int vflag) if (differentiation_flag == 1) gc->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); else gc->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); // extra per-atom energy/virial communication if (evflag_atom) { if (differentiation_flag == 1 && vflag_atom) gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); else if (differentiation_flag == 0) gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } // calculate the force on my particles @@ -817,10 +817,10 @@ void PPPM::allocate() // create ghost grid object for rho and electric field communication // also create 2 bufs for ghost grid cell comm, passed to GridComm methods - + gc = new GridComm(lmp,world,nx_pppm,ny_pppm,nz_pppm, - nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, - nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); + nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, + nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); gc->setup(ngc_buf1,ngc_buf2); @@ -988,7 +988,7 @@ void PPPM::set_grid_global() while (1) { // set grid dimensions - + nx_pppm = static_cast (xprd/h_x); ny_pppm = static_cast (yprd/h_y); nz_pppm = static_cast (zprd_slab/h_z); @@ -997,14 +997,14 @@ void PPPM::set_grid_global() if (ny_pppm <= 1) ny_pppm = 2; if (nz_pppm <= 1) nz_pppm = 2; - // estimate Kspace force error - + // estimate Kspace force error + double df_kspace = compute_df_kspace(); // break loop if the accuracy has been reached or // too many loops have been performed - count++; + count++; if (df_kspace <= accuracy) break; if (count > 500) error->all(FLERR, "Could not compute grid size"); @@ -1155,10 +1155,10 @@ double PPPM::compute_qopt() // loop over entire FFT grid // each proc calculates contributions from every Pth grid point - + bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm; int nxy_pppm = nx_pppm * ny_pppm; - + double qopt = 0.0; for (bigint i = me; i < ngridtotal; i += nprocs) { @@ -1174,46 +1174,46 @@ double PPPM::compute_qopt() if (sqk == 0.0) continue; sum1 = sum2 = sum3 = sum4 = 0.0; - + for (nx = -2; nx <= 2; nx++) { qx = unitkx*(kper+nx_pppm*nx); sx = exp(-0.25*square(qx/g_ewald)); argx = 0.5*qx*xprd/nx_pppm; wx = powsinxx(argx,twoorder); qx *= qx; - + for (ny = -2; ny <= 2; ny++) { - qy = unitky*(lper+ny_pppm*ny); - sy = exp(-0.25*square(qy/g_ewald)); - argy = 0.5*qy*yprd/ny_pppm; - wy = powsinxx(argy,twoorder); - qy *= qy; - - for (nz = -2; nz <= 2; nz++) { - qz = unitkz*(mper+nz_pppm*nz); - sz = exp(-0.25*square(qz/g_ewald)); - argz = 0.5*qz*zprd_slab/nz_pppm; - wz = powsinxx(argz,twoorder); - qz *= qz; - - dot2 = qx+qy+qz; - u1 = sx*sy*sz; - u2 = wx*wy*wz; - - sum1 += u1*u1/dot2*MY_4PI*MY_4PI; - sum2 += u1 * u2 * MY_4PI; - sum3 += u2; - sum4 += dot2*u2; - } + qy = unitky*(lper+ny_pppm*ny); + sy = exp(-0.25*square(qy/g_ewald)); + argy = 0.5*qy*yprd/ny_pppm; + wy = powsinxx(argy,twoorder); + qy *= qy; + + for (nz = -2; nz <= 2; nz++) { + qz = unitkz*(mper+nz_pppm*nz); + sz = exp(-0.25*square(qz/g_ewald)); + argz = 0.5*qz*zprd_slab/nz_pppm; + wz = powsinxx(argz,twoorder); + qz *= qz; + + dot2 = qx+qy+qz; + u1 = sx*sy*sz; + u2 = wx*wy*wz; + + sum1 += u1*u1/dot2*MY_4PI*MY_4PI; + sum2 += u1 * u2 * MY_4PI; + sum3 += u2; + sum4 += dot2*u2; + } } } - + sum2 *= sum2; qopt += sum1 - sum2/(sum3*sum4); } // sum qopt over all procs - + double qopt_all; MPI_Allreduce(&qopt,&qopt_all,1,MPI_DOUBLE,MPI_SUM,world); return qopt_all; @@ -1327,7 +1327,7 @@ void PPPM::set_grid_local() // global PPPM grid that I own without ghost cells // for slab PPPM, assign z grid as if it were not extended // both non-tiled and tiled proc layouts use 0-1 fractional sumdomain info - + if (comm->layout != Comm::LAYOUT_TILED) { nxlo_in = static_cast (comm->xsplit[comm->myloc[0]] * nx_pppm); nxhi_in = static_cast (comm->xsplit[comm->myloc[0]+1] * nx_pppm) - 1; @@ -1339,7 +1339,7 @@ void PPPM::set_grid_local() (comm->zsplit[comm->myloc[2]] * nz_pppm/slab_volfactor); nzhi_in = static_cast (comm->zsplit[comm->myloc[2]+1] * nz_pppm/slab_volfactor) - 1; - + } else { nxlo_in = static_cast (comm->mysplit[0][0] * nx_pppm); nxhi_in = static_cast (comm->mysplit[0][1] * nx_pppm) - 1; @@ -1437,7 +1437,7 @@ void PPPM::set_grid_local() // this is accomplished by nzhi_in = nzhi_out on +z end (no ghost cells) // also insure no other procs use ghost cells beyond +z limit // differnet logic for non-tiled vs tiled decomposition - + if (slabflag == 1) { if (comm->layout != Comm::LAYOUT_TILED) { if (comm->myloc[2] == comm->procgrid[2]-1) nzhi_in = nzhi_out = nz_pppm - 1; @@ -2634,7 +2634,7 @@ void PPPM::fieldforce_peratom() void PPPM::pack_forward_grid(int flag, void *vbuf, int nlist, int *list) { FFT_SCALAR *buf = (FFT_SCALAR *) vbuf; - + int n = 0; if (flag == FORWARD_IK) { @@ -2754,7 +2754,7 @@ void PPPM::unpack_forward_grid(int flag, void *vbuf, int nlist, int *list) void PPPM::pack_reverse_grid(int flag, void *vbuf, int nlist, int *list) { FFT_SCALAR *buf = (FFT_SCALAR *) vbuf; - + if (flag == REVERSE_RHO) { FFT_SCALAR *src = &density_brick[nzlo_out][nylo_out][nxlo_out]; for (int i = 0; i < nlist; i++) @@ -3083,7 +3083,7 @@ double PPPM::memory_usage() // two GridComm bufs bytes += (ngc_buf1 + ngc_buf2) * npergrid * sizeof(FFT_SCALAR); - + return bytes; } @@ -3141,7 +3141,7 @@ void PPPM::compute_group_group(int groupbit_A, int groupbit_B, int AA_flag) density_fft = density_A_fft; gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft(); // group B @@ -3150,7 +3150,7 @@ void PPPM::compute_group_group(int groupbit_A, int groupbit_B, int AA_flag) density_fft = density_B_fft; gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft(); // switch back pointers diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h index b7416a0a9c..7451af47b4 100644 --- a/src/KSPACE/pppm.h +++ b/src/KSPACE/pppm.h @@ -101,7 +101,7 @@ class PPPM : public KSpace { class FFT3d *fft1,*fft2; class Remap *remap; class GridComm *gc; - + FFT_SCALAR *gc_buf1,*gc_buf2; int ngc_buf1,ngc_buf2,npergrid; diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp index 081113ea0d..392d19336a 100644 --- a/src/KSPACE/pppm_cg.cpp +++ b/src/KSPACE/pppm_cg.cpp @@ -91,7 +91,7 @@ void PPPMCG::compute(int eflag, int vflag) ev_init(eflag,vflag); if (evflag_atom && !peratom_allocate_flag) allocate_peratom(); - + // if atom count has changed, update qsum and qsqsum if (atom->natoms != natoms_original) { @@ -158,7 +158,7 @@ void PPPMCG::compute(int eflag, int vflag) } // only need to rebuild this list after a neighbor list update - + if (neighbor->ago == 0) { num_charged = 0; for (int i = 0; i < atom->nlocal; ++i) { @@ -180,7 +180,7 @@ void PPPMCG::compute(int eflag, int vflag) // remap from 3d decomposition to FFT decomposition gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft(); // compute potential gradient on my FFT grid and @@ -195,20 +195,20 @@ void PPPMCG::compute(int eflag, int vflag) if (differentiation_flag == 1) gc->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); else gc->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); // extra per-atom energy/virial communication if (evflag_atom) { if (differentiation_flag == 1 && vflag_atom) gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); else if (differentiation_flag == 0) gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } // calculate the force on my particles diff --git a/src/KSPACE/pppm_dipole.cpp b/src/KSPACE/pppm_dipole.cpp index 312da6c304..1a87c81bee 100644 --- a/src/KSPACE/pppm_dipole.cpp +++ b/src/KSPACE/pppm_dipole.cpp @@ -444,9 +444,9 @@ void PPPMDipole::compute(int eflag, int vflag) // all procs communicate density values from their ghost cells // to fully sum contribution in their 3d bricks // remap from 3d decomposition to FFT decomposition - + gc_dipole->reverse_comm_kspace(this,3,sizeof(FFT_SCALAR),REVERSE_MU, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft_dipole(); // compute potential gradient on my FFT grid and @@ -460,13 +460,13 @@ void PPPMDipole::compute(int eflag, int vflag) // to fill ghost cells surrounding their 3d bricks gc_dipole->forward_comm_kspace(this,9,sizeof(FFT_SCALAR),FORWARD_MU, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); // extra per-atom energy/virial communication if (evflag_atom) gc->forward_comm_kspace(this,18,sizeof(FFT_SCALAR),FORWARD_MU_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); // calculate the force on my particles @@ -510,7 +510,7 @@ void PPPMDipole::compute(int eflag, int vflag) for (i = 0; i < nlocal; i++) { eatom[i] *= 0.5; eatom[i] -= (mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + - mu[i][2]*mu[i][2])*2.0*g3/3.0/MY_PIS; + mu[i][2]*mu[i][2])*2.0*g3/3.0/MY_PIS; eatom[i] *= qscale; } } @@ -2483,7 +2483,7 @@ int PPPMDipole::timing_3d(int n, double &time3d) double PPPMDipole::memory_usage() { double bytes = nmax*3 * sizeof(double); - + int nbrick = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) * (nzhi_out-nzlo_out+1); bytes += 6 * nfft_both * sizeof(double); // vg diff --git a/src/KSPACE/pppm_dipole.h b/src/KSPACE/pppm_dipole.h index f7a8b63930..eb49361842 100644 --- a/src/KSPACE/pppm_dipole.h +++ b/src/KSPACE/pppm_dipole.h @@ -74,7 +74,7 @@ class PPPMDipole : public PPPM { int only_dipole_flag; double musum,musqsum,mu2; - + double find_gewald_dipole(double, double, bigint, double, double); double newton_raphson_f_dipole(double, double, bigint, double, double); double derivf_dipole(double, double, bigint, double, double); diff --git a/src/KSPACE/pppm_dipole_spin.cpp b/src/KSPACE/pppm_dipole_spin.cpp index c8cebdfeef..4e33c11793 100644 --- a/src/KSPACE/pppm_dipole_spin.cpp +++ b/src/KSPACE/pppm_dipole_spin.cpp @@ -302,7 +302,7 @@ void PPPMDipoleSpin::compute(int eflag, int vflag) // remap from 3d decomposition to FFT decomposition gc_dipole->reverse_comm_kspace(this,3,sizeof(FFT_SCALAR),REVERSE_MU, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft_dipole(); // compute potential gradient on my FFT grid and @@ -316,13 +316,13 @@ void PPPMDipoleSpin::compute(int eflag, int vflag) // to fill ghost cells surrounding their 3d bricks gc_dipole->forward_comm_kspace(this,9,sizeof(FFT_SCALAR),FORWARD_MU, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); // extra per-atom energy/virial communication if (evflag_atom) gc->forward_comm_kspace(this,18,sizeof(FFT_SCALAR),FORWARD_MU_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); // calculate the force on my particles diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index 313f41d2cc..8f141bb695 100644 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -410,7 +410,7 @@ void PPPMDisp::init() int iteration = 0; if (function[0]) { - + GridComm *gctmp = NULL; while (order >= minorder) { @@ -455,7 +455,7 @@ void PPPMDisp::init() error->all(FLERR,"Coulomb PPPMDisp order has been reduced below minorder"); if (!overlap_allowed && !gctmp->ghost_adjacent()) error->all(FLERR,"PPPMDisp grid stencil extends " - "beyond nearest neighbor processor"); + "beyond nearest neighbor processor"); if (gctmp) delete gctmp; // adjust g_ewald @@ -491,7 +491,7 @@ void PPPMDisp::init() iteration = 0; if (function[1] + function[2] + function[3]) { - + GridComm *gctmp = NULL; while (order_6 >= minorder) { @@ -523,7 +523,7 @@ void PPPMDisp::init() nzlo_in_6,nzhi_in_6, nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6, nzlo_out_6,nzhi_out_6); - + int tmp1,tmp2; gctmp->setup(tmp1,tmp2); if (gctmp->ghost_adjacent()) break; @@ -537,7 +537,7 @@ void PPPMDisp::init() "reduced below minorder"); if (!overlap_allowed && !gctmp->ghost_adjacent()) error->all(FLERR,"Dispersion PPPMDisp grid stencil extends " - "beyond nearest neighbor processor"); + "beyond nearest neighbor processor"); if (gctmp) delete gctmp; // adjust g_ewald_6 @@ -638,7 +638,7 @@ void PPPMDisp::setup() double unitkz = (2.0*MY_PI/zprd_slab); //compute the virial coefficients and green functions - + if (function[0]) { delxinv = nx_pppm/xprd; @@ -874,7 +874,7 @@ void PPPMDisp::setup_grid() void PPPMDisp::compute(int eflag, int vflag) { int i; - + // set energy/virial flags // invoke allocate_peratom() if needed for first time @@ -889,7 +889,7 @@ void PPPMDisp::compute(int eflag, int vflag) boxlo = domain->boxlo_lamda; domain->x2lamda(atom->nlocal); } - + // extend size of per-atom arrays if necessary if (atom->nmax > nmax) { @@ -918,12 +918,12 @@ void PPPMDisp::compute(int eflag, int vflag) // perform calculations for coulomb interactions only particle_map_c(delxinv, delyinv, delzinv, shift, part2grid, nupper, nlower, - nxlo_out, nylo_out, nzlo_out, nxhi_out, nyhi_out, nzhi_out); + nxlo_out, nylo_out, nzlo_out, nxhi_out, nyhi_out, nzhi_out); make_rho_c(); - + gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft(nxlo_in, nylo_in, nzlo_in, nxhi_in, nyhi_in, nzhi_in, density_brick, density_fft, work1,remap); @@ -938,13 +938,13 @@ void PPPMDisp::compute(int eflag, int vflag) u_brick, v0_brick, v1_brick, v2_brick, v3_brick, v4_brick, v5_brick); gc->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); fieldforce_c_ad(); if (vflag_atom) - gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } else { poisson_ik(work1, work2, density_fft, fft1, fft2, @@ -957,31 +957,31 @@ void PPPMDisp::compute(int eflag, int vflag) u_brick, v0_brick, v1_brick, v2_brick, v3_brick, v4_brick, v5_brick); gc->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); fieldforce_c_ik(); if (evflag_atom) - gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } - + if (evflag_atom) fieldforce_c_peratom(); } if (function[1]) { - + // perform calculations for geometric mixing - + particle_map(delxinv_6, delyinv_6, delzinv_6, shift_6, part2grid_6, - nupper_6, nlower_6, + nupper_6, nlower_6, nxlo_out_6, nylo_out_6, nzlo_out_6, - nxhi_out_6, nyhi_out_6, nzhi_out_6); - + nxhi_out_6, nyhi_out_6, nzhi_out_6); + make_rho_g(); gc6->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO_G, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft(nxlo_in_6, nylo_in_6, nzlo_in_6, nxhi_in_6, nyhi_in_6, nzhi_in_6, density_brick_g, density_fft_g, work1_6,remap_6); @@ -994,16 +994,16 @@ void PPPMDisp::compute(int eflag, int vflag) energy_6, greensfn_6, virial_6, vg_6, vg2_6, u_brick_g, v0_brick_g, v1_brick_g, v2_brick_g, - v3_brick_g, v4_brick_g, v5_brick_g); + v3_brick_g, v4_brick_g, v5_brick_g); gc6->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD_G, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); fieldforce_g_ad(); if (vflag_atom) - gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_G, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_G, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } else { poisson_ik(work1_6, work2_6, density_fft_g, fft1_6, fft2_6, @@ -1014,34 +1014,34 @@ void PPPMDisp::compute(int eflag, int vflag) fkx_6, fky_6, fkz_6,fkx2_6, fky2_6, fkz2_6, vdx_brick_g, vdy_brick_g, vdz_brick_g, virial_6, vg_6, vg2_6, u_brick_g, v0_brick_g, v1_brick_g, v2_brick_g, - v3_brick_g, v4_brick_g, v5_brick_g); + v3_brick_g, v4_brick_g, v5_brick_g); gc6->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK_G, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); fieldforce_g_ik(); if (evflag_atom) - gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_G, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_G, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } - + if (evflag_atom) fieldforce_g_peratom(); } if (function[2]) { - + // perform calculations for arithmetic mixing - + particle_map(delxinv_6, delyinv_6, delzinv_6, shift_6, part2grid_6, - nupper_6, nlower_6, + nupper_6, nlower_6, nxlo_out_6, nylo_out_6, nzlo_out_6, - nxhi_out_6, nyhi_out_6, nzhi_out_6); - + nxhi_out_6, nyhi_out_6, nzhi_out_6); + make_rho_a(); gc->reverse_comm_kspace(this,7,sizeof(FFT_SCALAR),REVERSE_RHO_A, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft_a(); @@ -1053,31 +1053,31 @@ void PPPMDisp::compute(int eflag, int vflag) energy_6, greensfn_6, virial_6, vg_6, vg2_6, u_brick_a3, v0_brick_a3, v1_brick_a3, v2_brick_a3, - v3_brick_a3, v4_brick_a3, v5_brick_a3); + v3_brick_a3, v4_brick_a3, v5_brick_a3); poisson_2s_ad(density_fft_a0, density_fft_a6, u_brick_a0, v0_brick_a0, v1_brick_a0, v2_brick_a0, - v3_brick_a0, v4_brick_a0, v5_brick_a0, + v3_brick_a0, v4_brick_a0, v5_brick_a0, u_brick_a6, v0_brick_a6, v1_brick_a6, v2_brick_a6, - v3_brick_a6, v4_brick_a6, v5_brick_a6); + v3_brick_a6, v4_brick_a6, v5_brick_a6); poisson_2s_ad(density_fft_a1, density_fft_a5, u_brick_a1, v0_brick_a1, v1_brick_a1, v2_brick_a1, - v3_brick_a1, v4_brick_a1, v5_brick_a1, + v3_brick_a1, v4_brick_a1, v5_brick_a1, u_brick_a5, v0_brick_a5, v1_brick_a5, v2_brick_a5, - v3_brick_a5, v4_brick_a5, v5_brick_a5); + v3_brick_a5, v4_brick_a5, v5_brick_a5); poisson_2s_ad(density_fft_a2, density_fft_a4, u_brick_a2, v0_brick_a2, v1_brick_a2, v2_brick_a2, - v3_brick_a2, v4_brick_a2, v5_brick_a2, + v3_brick_a2, v4_brick_a2, v5_brick_a2, u_brick_a4, v0_brick_a4, v1_brick_a4, v2_brick_a4, - v3_brick_a4, v4_brick_a4, v5_brick_a4); + v3_brick_a4, v4_brick_a4, v5_brick_a4); gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_AD_A, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); fieldforce_a_ad(); if (evflag_atom) - gc6->forward_comm_kspace(this,42,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_A, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,42,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_A, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } else { poisson_ik(work1_6, work2_6, density_fft_a3, fft1_6, fft2_6, @@ -1088,55 +1088,55 @@ void PPPMDisp::compute(int eflag, int vflag) fkx_6, fky_6, fkz_6,fkx2_6, fky2_6, fkz2_6, vdx_brick_a3, vdy_brick_a3, vdz_brick_a3, virial_6, vg_6, vg2_6, u_brick_a3, v0_brick_a3, v1_brick_a3, v2_brick_a3, - v3_brick_a3, v4_brick_a3, v5_brick_a3); + v3_brick_a3, v4_brick_a3, v5_brick_a3); poisson_2s_ik(density_fft_a0, density_fft_a6, vdx_brick_a0, vdy_brick_a0, vdz_brick_a0, vdx_brick_a6, vdy_brick_a6, vdz_brick_a6, u_brick_a0, v0_brick_a0, v1_brick_a0, v2_brick_a0, - v3_brick_a0, v4_brick_a0, v5_brick_a0, + v3_brick_a0, v4_brick_a0, v5_brick_a0, u_brick_a6, v0_brick_a6, v1_brick_a6, v2_brick_a6, - v3_brick_a6, v4_brick_a6, v5_brick_a6); + v3_brick_a6, v4_brick_a6, v5_brick_a6); poisson_2s_ik(density_fft_a1, density_fft_a5, vdx_brick_a1, vdy_brick_a1, vdz_brick_a1, vdx_brick_a5, vdy_brick_a5, vdz_brick_a5, u_brick_a1, v0_brick_a1, v1_brick_a1, v2_brick_a1, - v3_brick_a1, v4_brick_a1, v5_brick_a1, + v3_brick_a1, v4_brick_a1, v5_brick_a1, u_brick_a5, v0_brick_a5, v1_brick_a5, v2_brick_a5, - v3_brick_a5, v4_brick_a5, v5_brick_a5); + v3_brick_a5, v4_brick_a5, v5_brick_a5); poisson_2s_ik(density_fft_a2, density_fft_a4, vdx_brick_a2, vdy_brick_a2, vdz_brick_a2, vdx_brick_a4, vdy_brick_a4, vdz_brick_a4, u_brick_a2, v0_brick_a2, v1_brick_a2, v2_brick_a2, - v3_brick_a2, v4_brick_a2, v5_brick_a2, + v3_brick_a2, v4_brick_a2, v5_brick_a2, u_brick_a4, v0_brick_a4, v1_brick_a4, v2_brick_a4, - v3_brick_a4, v4_brick_a4, v5_brick_a4); + v3_brick_a4, v4_brick_a4, v5_brick_a4); gc6->forward_comm_kspace(this,18,sizeof(FFT_SCALAR),FORWARD_IK_A, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); fieldforce_a_ik(); if (evflag_atom) - gc6->forward_comm_kspace(this,49,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_A, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,49,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_A, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } - + if (evflag_atom) fieldforce_a_peratom(); } if (function[3]) { - + // perform calculations if no mixing rule applies - + particle_map(delxinv_6, delyinv_6, delzinv_6, shift_6, part2grid_6, - nupper_6, nlower_6, + nupper_6, nlower_6, nxlo_out_6, nylo_out_6, nzlo_out_6, - nxhi_out_6, nyhi_out_6, nzhi_out_6); + nxhi_out_6, nyhi_out_6, nzhi_out_6); make_rho_none(); gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO_NONE, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft_none(); @@ -1151,13 +1151,13 @@ void PPPMDisp::compute(int eflag, int vflag) } gc6->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD_NONE, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); fieldforce_none_ad(); if (vflag_atom) - gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_NONE, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_NONE, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } else { int n = 0; @@ -1172,15 +1172,15 @@ void PPPMDisp::compute(int eflag, int vflag) } gc6->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK_NONE, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); fieldforce_none_ik(); if (evflag_atom) - gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_NONE, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_NONE, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } - + if (evflag_atom) fieldforce_none_peratom(); } @@ -1194,7 +1194,7 @@ void PPPMDisp::compute(int eflag, int vflag) // sum energy across procs and add in volume-dependent term const double qscale = force->qqrd2e * scale; - + if (eflag_global) { double energy_all; MPI_Allreduce(&energy_1,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); @@ -1234,7 +1234,7 @@ void PPPMDisp::compute(int eflag, int vflag) // coulomb self energy correction for (i = 0; i < atom->nlocal; i++) { eatom[i] -= qscale*g_ewald*q[i]*q[i]/MY_PIS + - qscale*MY_PI2*q[i]*qsum / (g_ewald*g_ewald*volume); + qscale*MY_PI2*q[i]*qsum / (g_ewald*g_ewald*volume); } } if (function[1] + function[2] + function[3]) { @@ -1254,7 +1254,7 @@ void PPPMDisp::compute(int eflag, int vflag) for (i = 0; i < atom->nlocal; i++) { tmp = atom->type[i]; for (int n = 0; n < 3; n++) - vatom[i][n] -= MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumi[tmp]; + vatom[i][n] -= MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumi[tmp]; } } } @@ -1665,7 +1665,7 @@ void _noopt PPPMDisp::allocate() memory->create2d_offset(rho_coeff,order,(1-order)/2,order/2,"pppm/disp:rho_coeff"); memory->create2d_offset(drho1d,3,-order/2,order/2,"pppm/disp:rho1d"); memory->create2d_offset(drho_coeff,order,(1-order)/2,order/2, - "pppm/disp:drho_coeff"); + "pppm/disp:drho_coeff"); memory->create(greensfn,nfft_both,"pppm/disp:greensfn"); memory->create(vg,nfft_both,6,"pppm/disp:vg"); @@ -1721,7 +1721,7 @@ void _noopt PPPMDisp::allocate() if (differentiation_flag) npergrid = 1; else npergrid = 3; - + memory->create(gc_buf1,npergrid*ngc_buf1,"pppm:gc_buf1"); memory->create(gc_buf2,npergrid*ngc_buf2,"pppm:gc_buf2"); } @@ -1741,10 +1741,10 @@ void _noopt PPPMDisp::allocate() memory->create(gf_b_6,order_6,"pppm/disp:gf_b_6"); memory->create2d_offset(rho1d_6,3,-order_6/2,order_6/2,"pppm/disp:rho1d_6"); memory->create2d_offset(rho_coeff_6,order_6,(1-order_6)/2,order_6/2, - "pppm/disp:rho_coeff_6"); + "pppm/disp:rho_coeff_6"); memory->create2d_offset(drho1d_6,3,-order_6/2,order_6/2,"pppm/disp:drho1d_6"); memory->create2d_offset(drho_coeff_6,order_6,(1-order_6)/2,order_6/2, - "pppm/disp:drho_coeff_6"); + "pppm/disp:drho_coeff_6"); memory->create(greensfn_6,nfft_both_6,"pppm/disp:greensfn_6"); memory->create(vg_6,nfft_both_6,6,"pppm/disp:vg_6"); @@ -1777,35 +1777,35 @@ void _noopt PPPMDisp::allocate() fft1_6 = new FFT3d(lmp,world,nx_pppm_6,ny_pppm_6,nz_pppm_6, - nxlo_fft_6,nxhi_fft_6,nylo_fft_6,nyhi_fft_6,nzlo_fft_6,nzhi_fft_6, - nxlo_fft_6,nxhi_fft_6,nylo_fft_6,nyhi_fft_6,nzlo_fft_6,nzhi_fft_6, - 0,0,&tmp,collective_flag); + nxlo_fft_6,nxhi_fft_6,nylo_fft_6,nyhi_fft_6,nzlo_fft_6,nzhi_fft_6, + nxlo_fft_6,nxhi_fft_6,nylo_fft_6,nyhi_fft_6,nzlo_fft_6,nzhi_fft_6, + 0,0,&tmp,collective_flag); fft2_6 = new FFT3d(lmp,world,nx_pppm_6,ny_pppm_6,nz_pppm_6, - nxlo_fft_6,nxhi_fft_6,nylo_fft_6,nyhi_fft_6,nzlo_fft_6,nzhi_fft_6, - nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, - 0,0,&tmp,collective_flag); + nxlo_fft_6,nxhi_fft_6,nylo_fft_6,nyhi_fft_6,nzlo_fft_6,nzhi_fft_6, + nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, + 0,0,&tmp,collective_flag); remap_6 = new Remap(lmp,world, - nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, - nxlo_fft_6,nxhi_fft_6,nylo_fft_6,nyhi_fft_6,nzlo_fft_6,nzhi_fft_6, - 1,0,0,FFT_PRECISION,collective_flag); + nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, + nxlo_fft_6,nxhi_fft_6,nylo_fft_6,nyhi_fft_6,nzlo_fft_6,nzhi_fft_6, + 1,0,0,FFT_PRECISION,collective_flag); // create ghost grid object for rho and electric field communication // also create 2 bufs for ghost grid cell comm, passed to GridComm methods gc6 = new GridComm(lmp,world,nx_pppm_6,ny_pppm_6,nz_pppm_6, - nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, - nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6); + nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, + nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6); gc6->setup(ngc6_buf1,ngc6_buf2); if (differentiation_flag) npergrid6 = 1; else npergrid6 = 7; - + memory->create(gc6_buf1,npergrid6*ngc6_buf1,"pppm:gc_buf1"); memory->create(gc6_buf2,npergrid6*ngc6_buf2,"pppm:gc_buf2"); } @@ -1825,35 +1825,35 @@ void _noopt PPPMDisp::allocate() memory->create(gf_b_6,order_6,"pppm/disp:gf_b_6"); memory->create2d_offset(rho1d_6,3,-order_6/2,order_6/2,"pppm/disp:rho1d_6"); memory->create2d_offset(rho_coeff_6,order_6,(1-order_6)/2,order_6/2, - "pppm/disp:rho_coeff_6"); + "pppm/disp:rho_coeff_6"); memory->create2d_offset(drho1d_6,3,-order_6/2,order_6/2,"pppm/disp:drho1d_6"); memory->create2d_offset(drho_coeff_6,order_6,(1-order_6)/2,order_6/2, - "pppm/disp:drho_coeff_6"); + "pppm/disp:drho_coeff_6"); memory->create(greensfn_6,nfft_both_6,"pppm/disp:greensfn_6"); memory->create(vg_6,nfft_both_6,6,"pppm/disp:vg_6"); memory->create(vg2_6,nfft_both_6,3,"pppm/disp:vg2_6"); memory->create3d_offset(density_brick_a0,nzlo_out_6, - nzhi_out_6,nylo_out_6,nyhi_out_6, + nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:density_brick_a0"); memory->create3d_offset(density_brick_a1,nzlo_out_6, - nzhi_out_6,nylo_out_6,nyhi_out_6, + nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:density_brick_a1"); memory->create3d_offset(density_brick_a2,nzlo_out_6,nzhi_out_6, - nylo_out_6,nyhi_out_6, + nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:density_brick_a2"); memory->create3d_offset(density_brick_a3,nzlo_out_6,nzhi_out_6, - nylo_out_6,nyhi_out_6, + nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:density_brick_a3"); memory->create3d_offset(density_brick_a4,nzlo_out_6,nzhi_out_6, - nylo_out_6,nyhi_out_6, + nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:density_brick_a4"); memory->create3d_offset(density_brick_a5,nzlo_out_6,nzhi_out_6, - nylo_out_6,nyhi_out_6, + nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:density_brick_a5"); memory->create3d_offset(density_brick_a6,nzlo_out_6,nzhi_out_6, - nylo_out_6,nyhi_out_6, + nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:density_brick_a6"); memory->create(density_fft_a0,nfft_both_6,"pppm/disp:density_fft_a0"); @@ -1959,17 +1959,17 @@ void _noopt PPPMDisp::allocate() // create ghost grid object for rho and electric field communication // also create 2 bufs for ghost grid cell comm, passed to GridComm methods - + gc6 = new GridComm(lmp,world,nx_pppm_6,ny_pppm_6,nz_pppm_6, - nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, - nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6); + nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, + nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6); gc6->setup(ngc6_buf1,ngc6_buf2); if (differentiation_flag) npergrid6 = 7; else npergrid6 = 18; - + memory->create(gc6_buf1,npergrid6*ngc6_buf1,"pppm:gc_buf1"); memory->create(gc6_buf2,npergrid6*ngc6_buf2,"pppm:gc_buf2"); } @@ -1989,21 +1989,21 @@ void _noopt PPPMDisp::allocate() memory->create(gf_b_6,order_6,"pppm/disp:gf_b_6"); memory->create2d_offset(rho1d_6,3,-order_6/2,order_6/2,"pppm/disp:rho1d_6"); memory->create2d_offset(rho_coeff_6,order_6,(1-order_6)/2,order_6/2, - "pppm/disp:rho_coeff_6"); + "pppm/disp:rho_coeff_6"); memory->create2d_offset(drho1d_6,3,-order_6/2,order_6/2,"pppm/disp:drho1d_6"); memory->create2d_offset(drho_coeff_6,order_6,(1-order_6)/2,order_6/2, - "pppm/disp:drho_coeff_6"); + "pppm/disp:drho_coeff_6"); memory->create(greensfn_6,nfft_both_6,"pppm/disp:greensfn_6"); memory->create(vg_6,nfft_both_6,6,"pppm/disp:vg_6"); memory->create(vg2_6,nfft_both_6,3,"pppm/disp:vg2_6"); memory->create4d_offset(density_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:density_brick_none"); if ( differentiation_flag == 1) { memory->create4d_offset(u_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:u_brick_none"); memory->create(sf_precoeff1_6,nfft_both_6,"pppm/disp:sf_precoeff1_6"); @@ -2015,17 +2015,17 @@ void _noopt PPPMDisp::allocate() } else { memory->create4d_offset(vdx_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:vdx_brick_none"); memory->create4d_offset(vdy_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:vdy_brick_none"); memory->create4d_offset(vdz_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:vdz_brick_none"); } memory->create(density_fft_none,nsplit_alloc,nfft_both_6, - "pppm/disp:density_fft_none"); + "pppm/disp:density_fft_none"); int tmp; @@ -2046,17 +2046,17 @@ void _noopt PPPMDisp::allocate() // create ghost grid object for rho and electric field communication // also create 2 bufs for ghost grid cell comm, passed to GridComm methods - + gc6 = new GridComm(lmp,world,nx_pppm_6,ny_pppm_6,nz_pppm_6, - nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, - nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6); + nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, + nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6); gc6->setup(ngc6_buf1,ngc6_buf2); if (differentiation_flag) npergrid6 = 1; else npergrid6 = 3; - + memory->create(gc6_buf1,npergrid6*ngc6_buf1,"pppm:gc_buf1"); memory->create(gc6_buf2,npergrid6*ngc6_buf2,"pppm:gc_buf2"); } @@ -2253,26 +2253,26 @@ void PPPMDisp::allocate_peratom() if (function[3]) { if (differentiation_flag != 1) memory->create4d_offset(u_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:u_brick_none"); memory->create4d_offset(v0_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:v0_brick_none"); memory->create4d_offset(v1_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:v1_brick_none"); memory->create4d_offset(v2_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:v2_brick_none"); memory->create4d_offset(v3_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:v3_brick_none"); memory->create4d_offset(v4_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:v4_brick_none"); memory->create4d_offset(v5_brick_none,nsplit_alloc, - nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, + nzlo_out_6,nzhi_out_6,nylo_out_6,nyhi_out_6, nxlo_out_6,nxhi_out_6,"pppm/disp:v5_brick_none"); // use same GC ghost grid object for peratom grid communication @@ -2604,7 +2604,7 @@ void PPPMDisp::set_grid() 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); @@ -2623,7 +2623,7 @@ void PPPMDisp::set_grid() count++; if (dfkspace <= accuracy) break; - + if (count > 500) error->all(FLERR, "Could not compute grid size"); h *= 0.95; h_x = h_y = h_z = h; @@ -2642,15 +2642,15 @@ void PPPMDisp::set_grid() ------------------------------------------------------------------------- */ void PPPMDisp::set_fft_parameters(int& nx_p,int& ny_p,int& nz_p, - int& nxlo_f,int& nylo_f,int& nzlo_f, - int& nxhi_f,int& nyhi_f,int& nzhi_f, - int& nxlo_i,int& nylo_i,int& nzlo_i, - int& nxhi_i,int& nyhi_i,int& nzhi_i, - int& nxlo_o,int& nylo_o,int& nzlo_o, - int& nxhi_o,int& nyhi_o,int& nzhi_o, - int& nlow, int& nupp, - int& ng, int& nf, int& nfb, - double& sft,double& sftone, int& ord) + int& nxlo_f,int& nylo_f,int& nzlo_f, + int& nxhi_f,int& nyhi_f,int& nzhi_f, + int& nxlo_i,int& nylo_i,int& nzlo_i, + int& nxhi_i,int& nyhi_i,int& nzhi_i, + int& nxlo_o,int& nylo_o,int& nzlo_o, + int& nxhi_o,int& nyhi_o,int& nzhi_o, + int& nlow, int& nupp, + int& ng, int& nf, int& nfb, + double& sft,double& sftone, int& ord) { // global indices of PPPM grid range from 0 to N-1 // nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of @@ -2988,7 +2988,7 @@ double PPPMDisp::compute_qopt_ik() // loop over entire FFT grid // each proc calculates contributions from every Pth grid point - + bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm; int nxy_pppm = nx_pppm * ny_pppm; @@ -3005,45 +3005,45 @@ double PPPMDisp::compute_qopt_ik() sqk = pow(unitkx*kper,2.0) + pow(unitky*lper,2.0) + pow(unitkz*mper,2.0); if (sqk == 0.0) continue; - + sum1 = sum2 = sum3 = 0.0; - + for (nx = -nbx; nx <= nbx; nx++) { qx = unitkx*(kper+nx_pppm*nx); sx = exp(-0.25*pow(qx/g_ewald,2.0)); wx = 1.0; argx = 0.5*qx*xprd/nx_pppm; if (argx != 0.0) wx = pow(sin(argx)/argx,order); - + for (ny = -nby; ny <= nby; ny++) { - qy = unitky*(lper+ny_pppm*ny); - sy = exp(-0.25*pow(qy/g_ewald,2.0)); - wy = 1.0; - argy = 0.5*qy*yprd/ny_pppm; - if (argy != 0.0) wy = pow(sin(argy)/argy,order); - - for (nz = -nbz; nz <= nbz; nz++) { - qz = unitkz*(mper+nz_pppm*nz); - sz = exp(-0.25*pow(qz/g_ewald,2.0)); - wz = 1.0; - argz = 0.5*qz*zprd_slab/nz_pppm; - if (argz != 0.0) wz = pow(sin(argz)/argz,order); - - dot1 = unitkx*kper*qx + unitky*lper*qy + unitkz*mper*qz; - dot2 = qx*qx+qy*qy+qz*qz; - u2 = pow(wx*wy*wz,2.0); - sum1 += sx*sy*sz*sx*sy*sz/dot2*4.0*4.0*MY_PI*MY_PI; - sum2 += u2*sx*sy*sz*4.0*MY_PI/dot2*dot1; - sum3 += u2; - } + qy = unitky*(lper+ny_pppm*ny); + sy = exp(-0.25*pow(qy/g_ewald,2.0)); + wy = 1.0; + argy = 0.5*qy*yprd/ny_pppm; + if (argy != 0.0) wy = pow(sin(argy)/argy,order); + + for (nz = -nbz; nz <= nbz; nz++) { + qz = unitkz*(mper+nz_pppm*nz); + sz = exp(-0.25*pow(qz/g_ewald,2.0)); + wz = 1.0; + argz = 0.5*qz*zprd_slab/nz_pppm; + if (argz != 0.0) wz = pow(sin(argz)/argz,order); + + dot1 = unitkx*kper*qx + unitky*lper*qy + unitkz*mper*qz; + dot2 = qx*qx+qy*qy+qz*qz; + u2 = pow(wx*wy*wz,2.0); + sum1 += sx*sy*sz*sx*sy*sz/dot2*4.0*4.0*MY_PI*MY_PI; + sum2 += u2*sx*sy*sz*4.0*MY_PI/dot2*dot1; + sum3 += u2; + } } } - + sum2 *= sum2; sum3 *= sum3*sqk; qopt += sum1 -sum2/sum3; } - + return qopt; } @@ -3079,7 +3079,7 @@ double PPPMDisp::compute_qopt_ad() // loop over entire FFT grid // each proc calculates contributions from every Pth grid point - + bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm; int nxy_pppm = nx_pppm * ny_pppm; @@ -3105,31 +3105,31 @@ double PPPMDisp::compute_qopt_ad() wx = 1.0; argx = 0.5*qx*xprd/nx_pppm; if (argx != 0.0) wx = pow(sin(argx)/argx,order); - - for (ny = -nby; ny <= nby; ny++) { - qy = unitky*(lper+ny_pppm*ny); - sy = exp(-0.25*pow(qy/g_ewald,2.0)); - wy = 1.0; - argy = 0.5*qy*yprd/ny_pppm; - if (argy != 0.0) wy = pow(sin(argy)/argy,order); - - for (nz = -nbz; nz <= nbz; nz++) { - qz = unitkz*(mper+nz_pppm*nz); - sz = exp(-0.25*pow(qz/g_ewald,2.0)); - wz = 1.0; - argz = 0.5*qz*zprd_slab/nz_pppm; - if (argz != 0.0) wz = pow(sin(argz)/argz,order); - dot2 = qx*qx+qy*qy+qz*qz; - u2 = pow(wx*wy*wz,2.0); - sum1 += sx*sy*sz*sx*sy*sz/dot2*4.0*4.0*MY_PI*MY_PI; - sum2 += sx*sy*sz * u2*4.0*MY_PI; - sum3 += u2; - sum4 += dot2*u2; - } + for (ny = -nby; ny <= nby; ny++) { + qy = unitky*(lper+ny_pppm*ny); + sy = exp(-0.25*pow(qy/g_ewald,2.0)); + wy = 1.0; + argy = 0.5*qy*yprd/ny_pppm; + if (argy != 0.0) wy = pow(sin(argy)/argy,order); + + for (nz = -nbz; nz <= nbz; nz++) { + qz = unitkz*(mper+nz_pppm*nz); + sz = exp(-0.25*pow(qz/g_ewald,2.0)); + wz = 1.0; + argz = 0.5*qz*zprd_slab/nz_pppm; + if (argz != 0.0) wz = pow(sin(argz)/argz,order); + + dot2 = qx*qx+qy*qy+qz*qz; + u2 = pow(wx*wy*wz,2.0); + sum1 += sx*sy*sz*sx*sy*sz/dot2*4.0*4.0*MY_PI*MY_PI; + sum2 += sx*sy*sz * u2*4.0*MY_PI; + sum3 += u2; + sum4 += dot2*u2; + } } } - + sum2 *= sum2; qopt += sum1 - sum2/(sum3*sum4); } @@ -3173,7 +3173,7 @@ double PPPMDisp::compute_qopt_6_ik() // loop over entire FFT grid // each proc calculates contributions from every Pth grid point - + bigint ngridtotal = (bigint) nx_pppm_6 * ny_pppm_6 * nz_pppm_6; int nxy_pppm_6 = nx_pppm_6 * ny_pppm_6; @@ -3190,7 +3190,7 @@ double PPPMDisp::compute_qopt_6_ik() sqk = pow(unitkx*kper,2.0) + pow(unitky*lper,2.0) + pow(unitkz*mper,2.0); if (sqk == 0.0) continue; - + sum1 = sum2 = sum3 = 0.0; for (nx = -nbx; nx <= nbx; nx++) { @@ -3199,32 +3199,32 @@ double PPPMDisp::compute_qopt_6_ik() wx = 1.0; argx = 0.5*qx*xprd/nx_pppm_6; if (argx != 0.0) wx = pow(sin(argx)/argx,order_6); - - for (ny = -nby; ny <= nby; ny++) { - qy = unitky*(lper+ny_pppm_6*ny); - sy = exp(-qy*qy*inv2ew*inv2ew); - wy = 1.0; - argy = 0.5*qy*yprd/ny_pppm_6; - if (argy != 0.0) wy = pow(sin(argy)/argy,order_6); - - for (nz = -nbz; nz <= nbz; nz++) { - qz = unitkz*(mper+nz_pppm_6*nz); - sz = exp(-qz*qz*inv2ew*inv2ew); - wz = 1.0; - argz = 0.5*qz*zprd_slab/nz_pppm_6; - if (argz != 0.0) wz = pow(sin(argz)/argz,order_6); - dot1 = unitkx*kper*qx + unitky*lper*qy + unitkz*mper*qz; - dot2 = qx*qx+qy*qy+qz*qz; - rtdot2 = sqrt(dot2); - term = (1-2*dot2*inv2ew*inv2ew)*sx*sy*sz + - 2*dot2*rtdot2*inv2ew*inv2ew*inv2ew*rtpi*erfc(rtdot2*inv2ew); - term *= g_ewald_6*g_ewald_6*g_ewald_6; - u2 = pow(wx*wy*wz,2.0); - sum1 += term*term*MY_PI*MY_PI*MY_PI/9.0 * dot2; - sum2 += -u2*term*MY_PI*rtpi/3.0*dot1; - sum3 += u2; - } + for (ny = -nby; ny <= nby; ny++) { + qy = unitky*(lper+ny_pppm_6*ny); + sy = exp(-qy*qy*inv2ew*inv2ew); + wy = 1.0; + argy = 0.5*qy*yprd/ny_pppm_6; + if (argy != 0.0) wy = pow(sin(argy)/argy,order_6); + + for (nz = -nbz; nz <= nbz; nz++) { + qz = unitkz*(mper+nz_pppm_6*nz); + sz = exp(-qz*qz*inv2ew*inv2ew); + wz = 1.0; + argz = 0.5*qz*zprd_slab/nz_pppm_6; + if (argz != 0.0) wz = pow(sin(argz)/argz,order_6); + + dot1 = unitkx*kper*qx + unitky*lper*qy + unitkz*mper*qz; + dot2 = qx*qx+qy*qy+qz*qz; + rtdot2 = sqrt(dot2); + term = (1-2*dot2*inv2ew*inv2ew)*sx*sy*sz + + 2*dot2*rtdot2*inv2ew*inv2ew*inv2ew*rtpi*erfc(rtdot2*inv2ew); + term *= g_ewald_6*g_ewald_6*g_ewald_6; + u2 = pow(wx*wy*wz,2.0); + sum1 += term*term*MY_PI*MY_PI*MY_PI/9.0 * dot2; + sum2 += -u2*term*MY_PI*rtpi/3.0*dot1; + sum3 += u2; + } } } sum2 *= sum2; @@ -3271,7 +3271,7 @@ double PPPMDisp::compute_qopt_6_ad() // loop over entire FFT grid // each proc calculates contributions from every Pth grid point - + bigint ngridtotal = (bigint) nx_pppm_6 * ny_pppm_6 * nz_pppm_6; int nxy_pppm_6 = nx_pppm_6 * ny_pppm_6; @@ -3288,7 +3288,7 @@ double PPPMDisp::compute_qopt_6_ad() sqk = pow(unitkx*kper,2.0) + pow(unitky*lper,2.0) + pow(unitkz*mper,2.0); if (sqk == 0.0) continue; - + sum1 = sum2 = sum3 = sum4 = 0.0; for (nx = -nbx; nx <= nbx; nx++) { @@ -3299,30 +3299,30 @@ double PPPMDisp::compute_qopt_6_ad() if (argx != 0.0) wx = pow(sin(argx)/argx,order_6); for (ny = -nby; ny <= nby; ny++) { - qy = unitky*(lper+ny_pppm_6*ny); - sy = exp(-qy*qy*inv2ew*inv2ew); - wy = 1.0; - argy = 0.5*qy*yprd/ny_pppm_6; - if (argy != 0.0) wy = pow(sin(argy)/argy,order_6); + qy = unitky*(lper+ny_pppm_6*ny); + sy = exp(-qy*qy*inv2ew*inv2ew); + wy = 1.0; + argy = 0.5*qy*yprd/ny_pppm_6; + if (argy != 0.0) wy = pow(sin(argy)/argy,order_6); - for (nz = -nbz; nz <= nbz; nz++) { - qz = unitkz*(mper+nz_pppm_6*nz); - sz = exp(-qz*qz*inv2ew*inv2ew); - wz = 1.0; - argz = 0.5*qz*zprd_slab/nz_pppm_6; - if (argz != 0.0) wz = pow(sin(argz)/argz,order_6); + for (nz = -nbz; nz <= nbz; nz++) { + qz = unitkz*(mper+nz_pppm_6*nz); + sz = exp(-qz*qz*inv2ew*inv2ew); + wz = 1.0; + argz = 0.5*qz*zprd_slab/nz_pppm_6; + if (argz != 0.0) wz = pow(sin(argz)/argz,order_6); - dot2 = qx*qx+qy*qy+qz*qz; - rtdot2 = sqrt(dot2); - term = (1-2*dot2*inv2ew*inv2ew)*sx*sy*sz + - 2*dot2*rtdot2*inv2ew*inv2ew*inv2ew*rtpi*erfc(rtdot2*inv2ew); - term *= g_ewald_6*g_ewald_6*g_ewald_6; - u2 = pow(wx*wy*wz,2.0); - sum1 += term*term*MY_PI*MY_PI*MY_PI/9.0 * dot2; - sum2 += -term*MY_PI*rtpi/3.0 * u2 * dot2; - sum3 += u2; - sum4 += dot2*u2; - } + dot2 = qx*qx+qy*qy+qz*qz; + rtdot2 = sqrt(dot2); + term = (1-2*dot2*inv2ew*inv2ew)*sx*sy*sz + + 2*dot2*rtdot2*inv2ew*inv2ew*inv2ew*rtpi*erfc(rtdot2*inv2ew); + term *= g_ewald_6*g_ewald_6*g_ewald_6; + u2 = pow(wx*wy*wz,2.0); + sum1 += term*term*MY_PI*MY_PI*MY_PI/9.0 * dot2; + sum2 += -term*MY_PI*rtpi/3.0 * u2 * dot2; + sum3 += u2; + sum4 += dot2*u2; + } } } sum2 *= sum2; @@ -3340,7 +3340,7 @@ double PPPMDisp::compute_qopt_6_ad() void PPPMDisp::set_grid_6() { // calculate csum - + if (!csumflag) calc_csum(); if (!gewaldflag_6) set_init_g6(); if (!gridflag_6) set_n_pppm_6(); @@ -3416,9 +3416,9 @@ void PPPMDisp::calc_csum() MPI_Allreduce(neach,neach_all,ntypes+1,MPI_INT,MPI_SUM,world); // copmute csumij and csumi - + double d1, d2; - + if (function[1]) { for (i=1; i<=ntypes; i++) { for (j=1; j<=ntypes; j++) { @@ -3430,7 +3430,7 @@ void PPPMDisp::calc_csum() } } } - + if (function[2]) { for (i=1; i<=ntypes; i++) { for (j=1; j<=ntypes; j++) { @@ -3444,7 +3444,7 @@ void PPPMDisp::calc_csum() } } } - + if (function[3]) { for (i=1; i<=ntypes; i++) { for (j=1; j<=ntypes; j++) { @@ -3469,7 +3469,7 @@ void PPPMDisp::calc_csum() void PPPMDisp::adjust_gewald_6() { // use Newton solver to find g_ewald_6 - + double dx; // start loop @@ -3549,13 +3549,13 @@ void PPPMDisp::set_init_g6() // if df_real > 0, repeat divide g_ewald_6 by 2 until df_real < 0 // else, repeat multiply g_ewald_6 by 2 until df_real > 0 // perform bisection for the last two values of - + double df_real; double g_ewald_old; double gmin, gmax; // check if there is a user defined accuracy - + double acc_rspace = accuracy; if (accuracy_real_6 > 0) acc_rspace = accuracy_real_6; @@ -3619,7 +3619,7 @@ void PPPMDisp::set_n_pppm_6() if (accuracy_kspace_6 > 0.0) acc_kspace = accuracy_kspace_6; // initial value for the grid spacing - + h = h_x = h_y = h_z = 4.0/g_ewald_6; // decrease grid spacing until required precision is obtained @@ -3661,7 +3661,7 @@ void PPPMDisp::set_n_pppm_6() // break loop if the accuracy has been reached or // too many loops have been performed - + if (df_kspace <= acc_kspace) break; if (count > 500) error->all(FLERR, "Could not compute grid size for Dispersion"); h *= 0.95; @@ -3771,10 +3771,10 @@ void PPPMDisp::compute_gf() ------------------------------------------------------------------------- */ void PPPMDisp::compute_sf_precoeff(int nxp, int nyp, int nzp, int ord, - int nxlo_ft, int nylo_ft, int nzlo_ft, - int nxhi_ft, int nyhi_ft, int nzhi_ft, - double *sf_pre1, double *sf_pre2, double *sf_pre3, - double *sf_pre4, double *sf_pre5, double *sf_pre6) + int nxlo_ft, int nylo_ft, int nzlo_ft, + int nxhi_ft, int nyhi_ft, int nzhi_ft, + double *sf_pre1, double *sf_pre2, double *sf_pre3, + double *sf_pre4, double *sf_pre5, double *sf_pre6) { int i,k,l,m,n; double *prd; @@ -4523,22 +4523,22 @@ void PPPMDisp::make_rho_none() ------------------------------------------------------------------------- */ void PPPMDisp::poisson_ik(FFT_SCALAR* wk1, FFT_SCALAR* wk2, - FFT_SCALAR* dfft, LAMMPS_NS::FFT3d* ft1, - LAMMPS_NS::FFT3d* ft2, - int nx_p, int ny_p, int nz_p, int nft, - int nxlo_ft, int nylo_ft, int nzlo_ft, - int nxhi_ft, int nyhi_ft, int nzhi_ft, - int nxlo_i, int nylo_i, int nzlo_i, - int nxhi_i, int nyhi_i, int nzhi_i, - double& egy, double* gfn, - double* kx, double* ky, double* kz, - double* kx2, double* ky2, double* kz2, - FFT_SCALAR*** vx_brick, FFT_SCALAR*** vy_brick, - FFT_SCALAR*** vz_brick, - double* vir, double** vcoeff, double** vcoeff2, - FFT_SCALAR*** u_pa, FFT_SCALAR*** v0_pa, - FFT_SCALAR*** v1_pa, FFT_SCALAR*** v2_pa, - FFT_SCALAR*** v3_pa, FFT_SCALAR*** v4_pa, FFT_SCALAR*** v5_pa) + FFT_SCALAR* dfft, LAMMPS_NS::FFT3d* ft1, + LAMMPS_NS::FFT3d* ft2, + int nx_p, int ny_p, int nz_p, int nft, + int nxlo_ft, int nylo_ft, int nzlo_ft, + int nxhi_ft, int nyhi_ft, int nzhi_ft, + int nxlo_i, int nylo_i, int nzlo_i, + int nxhi_i, int nyhi_i, int nzhi_i, + double& egy, double* gfn, + double* kx, double* ky, double* kz, + double* kx2, double* ky2, double* kz2, + FFT_SCALAR*** vx_brick, FFT_SCALAR*** vy_brick, + FFT_SCALAR*** vz_brick, + double* vir, double** vcoeff, double** vcoeff2, + FFT_SCALAR*** u_pa, FFT_SCALAR*** v0_pa, + FFT_SCALAR*** v1_pa, FFT_SCALAR*** v2_pa, + FFT_SCALAR*** v3_pa, FFT_SCALAR*** v4_pa, FFT_SCALAR*** v5_pa) { int i,j,k,n; @@ -4669,18 +4669,18 @@ void PPPMDisp::poisson_ik(FFT_SCALAR* wk1, FFT_SCALAR* wk2, ------------------------------------------------------------------------- */ void PPPMDisp::poisson_ad(FFT_SCALAR* wk1, FFT_SCALAR* wk2, - FFT_SCALAR* dfft, LAMMPS_NS::FFT3d* ft1,LAMMPS_NS::FFT3d* ft2, - int nx_p, int ny_p, int nz_p, int nft, - int nxlo_ft, int nylo_ft, int nzlo_ft, - int nxhi_ft, int nyhi_ft, int nzhi_ft, - int nxlo_i, int nylo_i, int nzlo_i, - int nxhi_i, int nyhi_i, int nzhi_i, - double& egy, double* gfn, - double* vir, double** vcoeff, double** vcoeff2, - FFT_SCALAR*** u_pa, FFT_SCALAR*** v0_pa, - FFT_SCALAR*** v1_pa, FFT_SCALAR*** v2_pa, - FFT_SCALAR*** v3_pa, FFT_SCALAR*** v4_pa, - FFT_SCALAR*** v5_pa) + FFT_SCALAR* dfft, LAMMPS_NS::FFT3d* ft1,LAMMPS_NS::FFT3d* ft2, + int nx_p, int ny_p, int nz_p, int nft, + int nxlo_ft, int nylo_ft, int nzlo_ft, + int nxhi_ft, int nyhi_ft, int nzhi_ft, + int nxlo_i, int nylo_i, int nzlo_i, + int nxhi_i, int nyhi_i, int nzhi_i, + double& egy, double* gfn, + double* vir, double** vcoeff, double** vcoeff2, + FFT_SCALAR*** u_pa, FFT_SCALAR*** v0_pa, + FFT_SCALAR*** v1_pa, FFT_SCALAR*** v2_pa, + FFT_SCALAR*** v3_pa, FFT_SCALAR*** v4_pa, + FFT_SCALAR*** v5_pa) { int i,j,k,n; double eng; @@ -4760,13 +4760,13 @@ void PPPMDisp::poisson_ad(FFT_SCALAR* wk1, FFT_SCALAR* wk2, ------------------------------------------------------------------------- */ void PPPMDisp:: poisson_peratom(FFT_SCALAR* wk1, FFT_SCALAR* wk2, LAMMPS_NS::FFT3d* ft2, - double** vcoeff, double** vcoeff2, int nft, - int nxlo_i, int nylo_i, int nzlo_i, - int nxhi_i, int nyhi_i, int nzhi_i, - FFT_SCALAR*** v0_pa, FFT_SCALAR*** v1_pa, - FFT_SCALAR*** v2_pa, - FFT_SCALAR*** v3_pa, FFT_SCALAR*** v4_pa, - FFT_SCALAR*** v5_pa) + double** vcoeff, double** vcoeff2, int nft, + int nxlo_i, int nylo_i, int nzlo_i, + int nxhi_i, int nyhi_i, int nzhi_i, + FFT_SCALAR*** v0_pa, FFT_SCALAR*** v1_pa, + FFT_SCALAR*** v2_pa, + FFT_SCALAR*** v3_pa, FFT_SCALAR*** v4_pa, + FFT_SCALAR*** v5_pa) { //v0 & v1 term int n, i, j, k; @@ -4834,14 +4834,14 @@ void PPPMDisp:: poisson_peratom(FFT_SCALAR* wk1, FFT_SCALAR* wk2, LAMMPS_NS::FFT void PPPMDisp:: poisson_2s_ik(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, - FFT_SCALAR*** vxbrick_1, FFT_SCALAR*** vybrick_1, FFT_SCALAR*** vzbrick_1, - FFT_SCALAR*** vxbrick_2, FFT_SCALAR*** vybrick_2, FFT_SCALAR*** vzbrick_2, - FFT_SCALAR*** u_pa_1, FFT_SCALAR*** v0_pa_1, - FFT_SCALAR*** v1_pa_1, FFT_SCALAR*** v2_pa_1, - FFT_SCALAR*** v3_pa_1, FFT_SCALAR*** v4_pa_1, FFT_SCALAR*** v5_pa_1, - FFT_SCALAR*** u_pa_2, FFT_SCALAR*** v0_pa_2, - FFT_SCALAR*** v1_pa_2, FFT_SCALAR*** v2_pa_2, - FFT_SCALAR*** v3_pa_2, FFT_SCALAR*** v4_pa_2, FFT_SCALAR*** v5_pa_2) + FFT_SCALAR*** vxbrick_1, FFT_SCALAR*** vybrick_1, FFT_SCALAR*** vzbrick_1, + FFT_SCALAR*** vxbrick_2, FFT_SCALAR*** vybrick_2, FFT_SCALAR*** vzbrick_2, + FFT_SCALAR*** u_pa_1, FFT_SCALAR*** v0_pa_1, + FFT_SCALAR*** v1_pa_1, FFT_SCALAR*** v2_pa_1, + FFT_SCALAR*** v3_pa_1, FFT_SCALAR*** v4_pa_1, FFT_SCALAR*** v5_pa_1, + FFT_SCALAR*** u_pa_2, FFT_SCALAR*** v0_pa_2, + FFT_SCALAR*** v1_pa_2, FFT_SCALAR*** v2_pa_2, + FFT_SCALAR*** v3_pa_2, FFT_SCALAR*** v4_pa_2, FFT_SCALAR*** v5_pa_2) { int i,j,k,n; @@ -4881,7 +4881,7 @@ poisson_2s_ik(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, n = 0; for (i = 0; i < nfft_6; i++) { eng = 2 * s2 * greensfn_6[i] * - (work1_6[n]*work2_6[n+1] - work1_6[n+1]*work2_6[n]); + (work1_6[n]*work2_6[n+1] - work1_6[n+1]*work2_6[n]); for (j = 0; j < 6; j++) virial_6[j] += eng*vg_6[i][j]; if (eflag_global)energy_6 += eng; n += 2; @@ -4894,9 +4894,9 @@ poisson_2s_ik(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, n += 2; } } - + // unify the two transformed vectors for efficient calculations later - + for ( i = 0; i < 2*nfft_6; i++) { work1_6[i] += work2_6[i]; } @@ -4997,9 +4997,9 @@ poisson_2s_ik(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, } if (vflag_atom) poisson_2s_peratom(v0_pa_1, v1_pa_1, v2_pa_1, - v3_pa_1, v4_pa_1, v5_pa_1, + v3_pa_1, v4_pa_1, v5_pa_1, v0_pa_2, v1_pa_2, v2_pa_2, - v3_pa_2, v4_pa_2, v5_pa_2); + v3_pa_2, v4_pa_2, v5_pa_2); } @@ -5010,13 +5010,13 @@ poisson_2s_ik(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, void PPPMDisp:: poisson_none_ik(int n1, int n2,FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, - FFT_SCALAR*** vxbrick_1, FFT_SCALAR*** vybrick_1, - FFT_SCALAR*** vzbrick_1, - FFT_SCALAR*** vxbrick_2, FFT_SCALAR*** vybrick_2, - FFT_SCALAR*** vzbrick_2, - FFT_SCALAR**** u_pa, FFT_SCALAR**** v0_pa, - FFT_SCALAR**** v1_pa, FFT_SCALAR**** v2_pa, - FFT_SCALAR**** v3_pa, FFT_SCALAR**** v4_pa, FFT_SCALAR**** v5_pa) + FFT_SCALAR*** vxbrick_1, FFT_SCALAR*** vybrick_1, + FFT_SCALAR*** vzbrick_1, + FFT_SCALAR*** vxbrick_2, FFT_SCALAR*** vybrick_2, + FFT_SCALAR*** vzbrick_2, + FFT_SCALAR**** u_pa, FFT_SCALAR**** v0_pa, + FFT_SCALAR**** v1_pa, FFT_SCALAR**** v2_pa, + FFT_SCALAR**** v3_pa, FFT_SCALAR**** v4_pa, FFT_SCALAR**** v5_pa) { int i,j,k,n; double eng; @@ -5025,7 +5025,7 @@ poisson_none_ik(int n1, int n2,FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, // transform charge/dispersion density (r -> k) // only one tansform required when energies and pressures not needed - + if (eflag_global + vflag_global == 0) { n = 0; for (i = 0; i < nfft_6; i++) { @@ -5037,7 +5037,7 @@ poisson_none_ik(int n1, int n2,FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, } // two transforms are required when energies and pressures are calculated - + else { n = 0; for (i = 0; i < nfft_6; i++) { @@ -5057,8 +5057,8 @@ poisson_none_ik(int n1, int n2,FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, n = 0; for (i = 0; i < nfft_6; i++) { eng = s2 * greensfn_6[i] * - (B[n1]*(work1_6[n]*work1_6[n] + work1_6[n+1]*work1_6[n+1]) + - B[n2]*(work2_6[n]*work2_6[n] + work2_6[n+1]*work2_6[n+1])); + (B[n1]*(work1_6[n]*work1_6[n] + work1_6[n+1]*work1_6[n+1]) + + B[n2]*(work2_6[n]*work2_6[n] + work2_6[n+1]*work2_6[n+1])); for (j = 0; j < 6; j++) virial_6[j] += eng*vg_6[i][j]; if (eflag_global)energy_6 += eng; n += 2; @@ -5068,14 +5068,14 @@ poisson_none_ik(int n1, int n2,FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, for (i = 0; i < nfft_6; i++) { energy_6 += s2 * greensfn_6[i] * - (B[n1]*(work1_6[n]*work1_6[n] + work1_6[n+1]*work1_6[n+1]) + - B[n2]*(work2_6[n]*work2_6[n] + work2_6[n+1]*work2_6[n+1])); + (B[n1]*(work1_6[n]*work1_6[n] + work1_6[n+1]*work1_6[n+1]) + + B[n2]*(work2_6[n]*work2_6[n] + work2_6[n+1]*work2_6[n+1])); n += 2; } } - + // unify the two transformed vectors for efficient calculations later - + for ( i = 0; i < 2*nfft_6; i++) { work1_6[i] += work2_6[i]; } @@ -5187,12 +5187,12 @@ poisson_none_ik(int n1, int n2,FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, void PPPMDisp:: poisson_2s_ad(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, - FFT_SCALAR*** u_pa_1, FFT_SCALAR*** v0_pa_1, - FFT_SCALAR*** v1_pa_1, FFT_SCALAR*** v2_pa_1, - FFT_SCALAR*** v3_pa_1, FFT_SCALAR*** v4_pa_1, FFT_SCALAR*** v5_pa_1, - FFT_SCALAR*** u_pa_2, FFT_SCALAR*** v0_pa_2, - FFT_SCALAR*** v1_pa_2, FFT_SCALAR*** v2_pa_2, - FFT_SCALAR*** v3_pa_2, FFT_SCALAR*** v4_pa_2, FFT_SCALAR*** v5_pa_2) + FFT_SCALAR*** u_pa_1, FFT_SCALAR*** v0_pa_1, + FFT_SCALAR*** v1_pa_1, FFT_SCALAR*** v2_pa_1, + FFT_SCALAR*** v3_pa_1, FFT_SCALAR*** v4_pa_1, FFT_SCALAR*** v5_pa_1, + FFT_SCALAR*** u_pa_2, FFT_SCALAR*** v0_pa_2, + FFT_SCALAR*** v1_pa_2, FFT_SCALAR*** v2_pa_2, + FFT_SCALAR*** v3_pa_2, FFT_SCALAR*** v4_pa_2, FFT_SCALAR*** v5_pa_2) { int i,j,k,n; double eng; @@ -5211,9 +5211,9 @@ poisson_2s_ad(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, fft1_6->compute(work1_6,work1_6,1); } - + // two transforms are required when energies and pressures are calculated - + else { n = 0; for (i = 0; i < nfft_6; i++) { @@ -5232,7 +5232,7 @@ poisson_2s_ad(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, n = 0; for (i = 0; i < nfft_6; i++) { eng = 2 * s2 * greensfn_6[i] * - (work1_6[n]*work2_6[n+1] - work1_6[n+1]*work2_6[n]); + (work1_6[n]*work2_6[n+1] - work1_6[n+1]*work2_6[n]); for (j = 0; j < 6; j++) virial_6[j] += eng*vg_6[i][j]; if (eflag_global)energy_6 += eng; n += 2; @@ -5277,9 +5277,9 @@ poisson_2s_ad(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, } if (vflag_atom) poisson_2s_peratom(v0_pa_1, v1_pa_1, v2_pa_1, - v3_pa_1, v4_pa_1, v5_pa_1, + v3_pa_1, v4_pa_1, v5_pa_1, v0_pa_2, v1_pa_2, v2_pa_2, - v3_pa_2, v4_pa_2, v5_pa_2); + v3_pa_2, v4_pa_2, v5_pa_2); } /* ---------------------------------------------------------------------- @@ -5289,9 +5289,9 @@ poisson_2s_ad(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, void PPPMDisp:: poisson_none_ad(int n1, int n2, FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, - FFT_SCALAR*** u_pa_1, FFT_SCALAR*** u_pa_2, - FFT_SCALAR**** v0_pa, FFT_SCALAR**** v1_pa, FFT_SCALAR**** v2_pa, - FFT_SCALAR**** v3_pa, FFT_SCALAR**** v4_pa, FFT_SCALAR**** v5_pa) + FFT_SCALAR*** u_pa_1, FFT_SCALAR*** u_pa_2, + FFT_SCALAR**** v0_pa, FFT_SCALAR**** v1_pa, FFT_SCALAR**** v2_pa, + FFT_SCALAR**** v3_pa, FFT_SCALAR**** v4_pa, FFT_SCALAR**** v5_pa) { int i,j,k,n; double eng; @@ -5310,9 +5310,9 @@ poisson_none_ad(int n1, int n2, FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, fft1_6->compute(work1_6,work1_6,1); } - + // two transforms are required when energies and pressures are calculated - + else { n = 0; for (i = 0; i < nfft_6; i++) { @@ -5331,8 +5331,8 @@ poisson_none_ad(int n1, int n2, FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, n = 0; for (i = 0; i < nfft_6; i++) { eng = s2 * greensfn_6[i] * - (B[n1]*(work1_6[n]*work1_6[n] + work1_6[n+1]*work1_6[n+1]) + - B[n2]*(work2_6[n]*work2_6[n] + work2_6[n+1]*work2_6[n+1])); + (B[n1]*(work1_6[n]*work1_6[n] + work1_6[n+1]*work1_6[n+1]) + + B[n2]*(work2_6[n]*work2_6[n] + work2_6[n+1]*work2_6[n+1])); for (j = 0; j < 6; j++) virial_6[j] += eng*vg_6[i][j]; if (eflag_global)energy_6 += eng; n += 2; @@ -5342,14 +5342,14 @@ poisson_none_ad(int n1, int n2, FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, for (i = 0; i < nfft_6; i++) { energy_6 += s2 * greensfn_6[i] * - (B[n1]*(work1_6[n]*work1_6[n] + work1_6[n+1]*work1_6[n+1]) + - B[n2]*(work2_6[n]*work2_6[n] + work2_6[n+1]*work2_6[n+1])); + (B[n1]*(work1_6[n]*work1_6[n] + work1_6[n+1]*work1_6[n+1]) + + B[n2]*(work2_6[n]*work2_6[n] + work2_6[n+1]*work2_6[n+1])); n += 2; } } - + // unify the two transformed vectors for efficient calculations later - + for ( i = 0; i < 2*nfft_6; i++) { work1_6[i] += work2_6[i]; } @@ -5391,12 +5391,12 @@ poisson_none_ad(int n1, int n2, FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2, void PPPMDisp:: poisson_2s_peratom(FFT_SCALAR*** v0_pa_1, FFT_SCALAR*** v1_pa_1, FFT_SCALAR*** v2_pa_1, - FFT_SCALAR*** v3_pa_1, FFT_SCALAR*** v4_pa_1, FFT_SCALAR*** v5_pa_1, - FFT_SCALAR*** v0_pa_2, FFT_SCALAR*** v1_pa_2, FFT_SCALAR*** v2_pa_2, - FFT_SCALAR*** v3_pa_2, FFT_SCALAR*** v4_pa_2, FFT_SCALAR*** v5_pa_2) + FFT_SCALAR*** v3_pa_1, FFT_SCALAR*** v4_pa_1, FFT_SCALAR*** v5_pa_1, + FFT_SCALAR*** v0_pa_2, FFT_SCALAR*** v1_pa_2, FFT_SCALAR*** v2_pa_2, + FFT_SCALAR*** v3_pa_2, FFT_SCALAR*** v4_pa_2, FFT_SCALAR*** v5_pa_2) { //Compute first virial term v0 - + int n, i, j, k; n = 0; @@ -5518,14 +5518,14 @@ poisson_2s_peratom(FFT_SCALAR*** v0_pa_1, FFT_SCALAR*** v1_pa_1, FFT_SCALAR*** v void PPPMDisp:: poisson_none_peratom(int n1, int n2, - FFT_SCALAR*** v0_pa_1, FFT_SCALAR*** v1_pa_1, - FFT_SCALAR*** v2_pa_1, - FFT_SCALAR*** v3_pa_1, FFT_SCALAR*** v4_pa_1, - FFT_SCALAR*** v5_pa_1, - FFT_SCALAR*** v0_pa_2, FFT_SCALAR*** v1_pa_2, - FFT_SCALAR*** v2_pa_2, - FFT_SCALAR*** v3_pa_2, FFT_SCALAR*** v4_pa_2, - FFT_SCALAR*** v5_pa_2) + FFT_SCALAR*** v0_pa_1, FFT_SCALAR*** v1_pa_1, + FFT_SCALAR*** v2_pa_1, + FFT_SCALAR*** v3_pa_1, FFT_SCALAR*** v4_pa_1, + FFT_SCALAR*** v5_pa_1, + FFT_SCALAR*** v0_pa_2, FFT_SCALAR*** v1_pa_2, + FFT_SCALAR*** v2_pa_2, + FFT_SCALAR*** v3_pa_2, FFT_SCALAR*** v4_pa_2, + FFT_SCALAR*** v5_pa_2) { //Compute first virial term v0 int n, i, j, k; @@ -6339,7 +6339,7 @@ void PPPMDisp::fieldforce_a_ad() sf *= 4*lj0*lj6 + 4*lj1*lj5 + 4*lj2*lj4 + 2*lj3*lj3; if (slabflag != 2) f[i][2] += lj0*ekz0 + lj1*ekz1 + - lj2*ekz2 + lj3*ekz3 + lj4*ekz4 + lj5*ekz5 + lj6*ekz6 - sf; + lj2*ekz2 + lj3*ekz3 + lj4*ekz4 + lj5*ekz5 + lj6*ekz6 - sf; } } @@ -8248,7 +8248,7 @@ int PPPMDisp::timing_3d(int n, double &time3d) double PPPMDisp::memory_usage() { double bytes = nmax*3 * sizeof(double); - + int mixing = 1; int diff = 3; //depends on differentiation int per = 7; //depends on per atom calculations @@ -8273,7 +8273,7 @@ double PPPMDisp::memory_usage() int nbrick = (nxhi_out_6-nxlo_out_6+1) * (nyhi_out_6-nylo_out_6+1) * (nzhi_out_6-nzlo_out_6+1); // density_brick + vd_brick + per atom bricks - bytes += (1 + diff + per ) * nbrick * sizeof(FFT_SCALAR) * mixing; + bytes += (1 + diff + per ) * nbrick * sizeof(FFT_SCALAR) * mixing; bytes += 6 * nfft_both_6 * sizeof(double); // vg bytes += nfft_both_6 * sizeof(double); // greensfn // density_FFT, work1, work2 diff --git a/src/KSPACE/pppm_stagger.cpp b/src/KSPACE/pppm_stagger.cpp index 837644f0e3..f71529ae83 100644 --- a/src/KSPACE/pppm_stagger.cpp +++ b/src/KSPACE/pppm_stagger.cpp @@ -157,7 +157,7 @@ void PPPMStagger::compute(int eflag, int vflag) // remap from 3d decomposition to FFT decomposition gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft(); // compute potential gradient on my FFT grid and @@ -172,20 +172,20 @@ void PPPMStagger::compute(int eflag, int vflag) if (differentiation_flag == 1) gc->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); else gc->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); // extra per-atom energy/virial communication if (evflag_atom) { if (differentiation_flag == 1 && vflag_atom) - gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); else if (differentiation_flag == 0) - gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } // calculate the force on my particles @@ -299,10 +299,10 @@ double PPPMStagger::compute_qopt() // loop over entire FFT grid // each proc calculates contributions from every Pth grid point - + bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm; int nxy_pppm = nx_pppm * ny_pppm; - + double qopt = 0.0; for (bigint i = me; i < ngridtotal; i += nprocs) { @@ -336,31 +336,31 @@ double PPPMStagger::compute_qopt() wx = powsinxx(argx,twoorder); for (ny = -nby; ny <= nby; ny++) { - qy = unitky*(lper+ny_pppm*ny); - sy = exp(-0.25*square(qy/g_ewald)); - argy = 0.5*qy*yprd/ny_pppm; - wy = powsinxx(argy,twoorder); + qy = unitky*(lper+ny_pppm*ny); + sy = exp(-0.25*square(qy/g_ewald)); + argy = 0.5*qy*yprd/ny_pppm; + wy = powsinxx(argy,twoorder); - for (nz = -nbz; nz <= nbz; nz++) { - qz = unitkz*(mper+nz_pppm*nz); - sz = exp(-0.25*square(qz/g_ewald)); - argz = 0.5*qz*zprd_slab/nz_pppm; - wz = powsinxx(argz,twoorder); + for (nz = -nbz; nz <= nbz; nz++) { + qz = unitkz*(mper+nz_pppm*nz); + sz = exp(-0.25*square(qz/g_ewald)); + argz = 0.5*qz*zprd_slab/nz_pppm; + wz = powsinxx(argz,twoorder); - dot1 = unitkx*kper*qx + unitky*lper*qy + unitkz*mper*qz; - dot2 = qx*qx + qy*qy + qz*qz; - u1 = sx*sy*sz; - u2 = wx*wy*wz; - u3 = numerator*u1*u2*dot1; - sum1 += u1*u1*MY_4PI*MY_4PI/dot2; - sum2 += u3*u3/dot2; - } + dot1 = unitkx*kper*qx + unitky*lper*qy + unitkz*mper*qz; + dot2 = qx*qx + qy*qy + qz*qz; + u1 = sx*sy*sz; + u2 = wx*wy*wz; + u3 = numerator*u1*u2*dot1; + sum1 += u1*u1*MY_4PI*MY_4PI/dot2; + sum2 += u3*u3/dot2; + } } } - + qopt += sum1 - sum2/denominator; } - + double qopt_all; MPI_Allreduce(&qopt,&qopt_all,1,MPI_DOUBLE,MPI_SUM,world); return qopt_all; @@ -395,10 +395,10 @@ double PPPMStagger::compute_qopt_ad() // loop over entire FFT grid // each proc calculates contributions from every Pth grid point - + bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm; int nxy_pppm = nx_pppm * ny_pppm; - + double qopt = 0.0; for (bigint i = me; i < ngridtotal; i += nprocs) { @@ -414,7 +414,7 @@ double PPPMStagger::compute_qopt_ad() if (sqk == 0.0) continue; sum1 = sum2 = sum3 = sum4 = sum5 = sum6 = 0.0; - + for (nx = -nbx; nx <= nbx; nx++) { qx = unitkx*(kper+nx_pppm*nx); sx = exp(-0.25*square(qx/g_ewald)); @@ -422,30 +422,30 @@ double PPPMStagger::compute_qopt_ad() wx = powsinxx(argx,twoorder); for (ny = -nby; ny <= nby; ny++) { - qy = unitky*(lper+ny_pppm*ny); - sy = exp(-0.25*square(qy/g_ewald)); - argy = 0.5*qy*yprd/ny_pppm; - wy = powsinxx(argy,twoorder); - - for (nz = -nbz; nz <= nbz; nz++) { - qz = unitkz*(mper+nz_pppm*nz); - sz = exp(-0.25*square(qz/g_ewald)); - argz = 0.5*qz*zprd_slab/nz_pppm; - wz = powsinxx(argz,twoorder); - - dot2 = qx*qx + qy*qy + qz*qz; - u1 = sx*sy*sz; - u2 = wx*wy*wz; - sum1 += u1*u1/dot2*MY_4PI*MY_4PI; - sum2 += u1*u1*u2*u2*MY_4PI*MY_4PI; - sum3 += u2; - sum4 += dot2*u2; - sum5 += u2*powint(-1.0,nx+ny+nz); - sum6 += dot2*u2*powint(-1.0,nx+ny+nz); - } + qy = unitky*(lper+ny_pppm*ny); + sy = exp(-0.25*square(qy/g_ewald)); + argy = 0.5*qy*yprd/ny_pppm; + wy = powsinxx(argy,twoorder); + + for (nz = -nbz; nz <= nbz; nz++) { + qz = unitkz*(mper+nz_pppm*nz); + sz = exp(-0.25*square(qz/g_ewald)); + argz = 0.5*qz*zprd_slab/nz_pppm; + wz = powsinxx(argz,twoorder); + + dot2 = qx*qx + qy*qy + qz*qz; + u1 = sx*sy*sz; + u2 = wx*wy*wz; + sum1 += u1*u1/dot2*MY_4PI*MY_4PI; + sum2 += u1*u1*u2*u2*MY_4PI*MY_4PI; + sum3 += u2; + sum4 += dot2*u2; + sum5 += u2*powint(-1.0,nx+ny+nz); + sum6 += dot2*u2*powint(-1.0,nx+ny+nz); + } } } - + qopt += sum1 - sum2/(0.5*(sum3*sum4 + sum5*sum6)); } diff --git a/src/MISC/fix_gld.cpp b/src/MISC/fix_gld.cpp index ab601ae8cb..0e4c61813b 100644 --- a/src/MISC/fix_gld.cpp +++ b/src/MISC/fix_gld.cpp @@ -568,7 +568,7 @@ void FixGLD::unpack_restart(int nlocal, int nth) // skip to the nth set of extended variables // unpack the Nth first values this way because other fixes pack them - + int m = 0; for (int i = 0; i< nth; i++) m += static_cast (extra[nlocal][m]); m++; diff --git a/src/MISC/fix_ttm.cpp b/src/MISC/fix_ttm.cpp index 529914ec34..896eb24b7c 100644 --- a/src/MISC/fix_ttm.cpp +++ b/src/MISC/fix_ttm.cpp @@ -674,7 +674,7 @@ void FixTTM::unpack_restart(int nlocal, int nth) // skip to Nth set of extra values // unpack the Nth first values this way because other fixes pack them - + int m = 0; for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); m++; diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp index d3b071ebc6..5d95ae73c2 100644 --- a/src/MOLECULE/fix_cmap.cpp +++ b/src/MOLECULE/fix_cmap.cpp @@ -1307,7 +1307,7 @@ void FixCMAP::unpack_restart(int nlocal, int nth) // skip to Nth set of extra values // unpack the Nth first values this way because other fixes pack them - + int n = 0; for (int i = 0; i < nth; i++) n += static_cast (extra[nlocal][n]); diff --git a/src/USER-INTEL/pppm_disp_intel.cpp b/src/USER-INTEL/pppm_disp_intel.cpp index 97e5c57d6e..3229c462aa 100644 --- a/src/USER-INTEL/pppm_disp_intel.cpp +++ b/src/USER-INTEL/pppm_disp_intel.cpp @@ -173,9 +173,9 @@ void PPPMDispIntel::compute(int eflag, int vflag) return; } #endif - + int i; - + // convert atoms from box to lamda coords ev_init(eflag,vflag); @@ -292,7 +292,7 @@ void PPPMDispIntel::compute(int eflag, int vflag) } gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft(nxlo_in, nylo_in, nzlo_in, nxhi_in, nyhi_in, nzhi_in, density_brick, density_fft, work1,remap); @@ -306,7 +306,7 @@ void PPPMDispIntel::compute(int eflag, int vflag) v1_brick, v2_brick, v3_brick, v4_brick, v5_brick); gc->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); if (fix->precision() == FixIntel::PREC_MODE_MIXED) { fieldforce_c_ad(fix->get_mixed_buffers()); @@ -317,8 +317,8 @@ void PPPMDispIntel::compute(int eflag, int vflag) } if (vflag_atom) - gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } else { poisson_ik(work1, work2, density_fft, fft1, fft2, @@ -331,7 +331,7 @@ void PPPMDispIntel::compute(int eflag, int vflag) v5_brick); gc->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); if (fix->precision() == FixIntel::PREC_MODE_MIXED) { fieldforce_c_ik(fix->get_mixed_buffers()); @@ -342,14 +342,14 @@ void PPPMDispIntel::compute(int eflag, int vflag) } if (evflag_atom) - gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } if (evflag_atom) fieldforce_c_peratom(); } if (function[1]) { - + //perform calculations for geometric mixing if (fix->precision() == FixIntel::PREC_MODE_MIXED) { @@ -376,7 +376,7 @@ void PPPMDispIntel::compute(int eflag, int vflag) } gc6->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO_G, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft(nxlo_in_6, nylo_in_6, nzlo_in_6, nxhi_in_6, nyhi_in_6, nzhi_in_6, density_brick_g, density_fft_g, work1_6,remap_6); @@ -391,19 +391,19 @@ void PPPMDispIntel::compute(int eflag, int vflag) v2_brick_g, v3_brick_g, v4_brick_g, v5_brick_g); gc6->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD_G, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); if (fix->precision() == FixIntel::PREC_MODE_MIXED) { - fieldforce_g_ad(fix->get_mixed_buffers()); + fieldforce_g_ad(fix->get_mixed_buffers()); } else if (fix->precision() == FixIntel::PREC_MODE_DOUBLE) { - fieldforce_g_ad(fix->get_double_buffers()); + fieldforce_g_ad(fix->get_double_buffers()); } else { - fieldforce_g_ad(fix->get_single_buffers()); + fieldforce_g_ad(fix->get_single_buffers()); } if (vflag_atom) - gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_G, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_G, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } else { poisson_ik(work1_6, work2_6, density_fft_g, fft1_6, fft2_6, @@ -416,19 +416,19 @@ void PPPMDispIntel::compute(int eflag, int vflag) v1_brick_g, v2_brick_g, v3_brick_g, v4_brick_g, v5_brick_g); gc6->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK_G, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); if (fix->precision() == FixIntel::PREC_MODE_MIXED) { - fieldforce_g_ik(fix->get_mixed_buffers()); + fieldforce_g_ik(fix->get_mixed_buffers()); } else if (fix->precision() == FixIntel::PREC_MODE_DOUBLE) { - fieldforce_g_ik(fix->get_double_buffers()); + fieldforce_g_ik(fix->get_double_buffers()); } else { - fieldforce_g_ik(fix->get_single_buffers()); + fieldforce_g_ik(fix->get_single_buffers()); } if (evflag_atom) - gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_G, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_G, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } if (evflag_atom) fieldforce_g_peratom(); @@ -461,7 +461,7 @@ void PPPMDispIntel::compute(int eflag, int vflag) } gc->reverse_comm_kspace(this,7,sizeof(FFT_SCALAR),REVERSE_RHO_A, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft_a(); @@ -487,19 +487,19 @@ void PPPMDispIntel::compute(int eflag, int vflag) v2_brick_a4, v3_brick_a4, v4_brick_a4, v5_brick_a4); gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_AD_A, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); if (fix->precision() == FixIntel::PREC_MODE_MIXED) { - fieldforce_a_ad(fix->get_mixed_buffers()); + fieldforce_a_ad(fix->get_mixed_buffers()); } else if (fix->precision() == FixIntel::PREC_MODE_DOUBLE) { - fieldforce_a_ad(fix->get_double_buffers()); + fieldforce_a_ad(fix->get_double_buffers()); } else { - fieldforce_a_ad(fix->get_single_buffers()); + fieldforce_a_ad(fix->get_single_buffers()); } if (evflag_atom) - gc6->forward_comm_kspace(this,42,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_A, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,42,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_A, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } else { poisson_ik(work1_6, work2_6, density_fft_a3, fft1_6, fft2_6, @@ -530,7 +530,7 @@ void PPPMDispIntel::compute(int eflag, int vflag) v3_brick_a4, v4_brick_a4, v5_brick_a4); gc6->forward_comm_kspace(this,18,sizeof(FFT_SCALAR),FORWARD_IK_A, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); if (fix->precision() == FixIntel::PREC_MODE_MIXED) { fieldforce_a_ik(fix->get_mixed_buffers()); @@ -541,15 +541,15 @@ void PPPMDispIntel::compute(int eflag, int vflag) } if (evflag_atom) - gc6->forward_comm_kspace(this,49,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_A, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,49,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_A, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } if (evflag_atom) fieldforce_a_peratom(); } if (function[3]) { - + // perform calculations if no mixing rule applies if (fix->precision() == FixIntel::PREC_MODE_MIXED) { @@ -576,7 +576,7 @@ void PPPMDispIntel::compute(int eflag, int vflag) } gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO_NONE, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); brick2fft_none(); @@ -592,19 +592,19 @@ void PPPMDispIntel::compute(int eflag, int vflag) } gc6->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD_NONE, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); if (fix->precision() == FixIntel::PREC_MODE_MIXED) { - fieldforce_none_ad(fix->get_mixed_buffers()); + fieldforce_none_ad(fix->get_mixed_buffers()); } else if (fix->precision() == FixIntel::PREC_MODE_DOUBLE) { - fieldforce_none_ad(fix->get_double_buffers()); + fieldforce_none_ad(fix->get_double_buffers()); } else { - fieldforce_none_ad(fix->get_single_buffers()); + fieldforce_none_ad(fix->get_single_buffers()); } if (vflag_atom) - gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_NONE, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_NONE, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } else { int n = 0; @@ -621,19 +621,19 @@ void PPPMDispIntel::compute(int eflag, int vflag) } gc6->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK_NONE, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc_buf1,gc_buf2,MPI_FFT_SCALAR); if (fix->precision() == FixIntel::PREC_MODE_MIXED) { - fieldforce_none_ik(fix->get_mixed_buffers()); + fieldforce_none_ik(fix->get_mixed_buffers()); } else if (fix->precision() == FixIntel::PREC_MODE_DOUBLE) { - fieldforce_none_ik(fix->get_double_buffers()); + fieldforce_none_ik(fix->get_double_buffers()); } else { - fieldforce_none_ik(fix->get_single_buffers()); + fieldforce_none_ik(fix->get_single_buffers()); } if (evflag_atom) - gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_NONE, - gc_buf1,gc_buf2,MPI_FFT_SCALAR); + gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_NONE, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); } if (evflag_atom) fieldforce_none_peratom(); diff --git a/src/USER-MISC/fix_gle.cpp b/src/USER-MISC/fix_gle.cpp index 8459ddf29b..87dbb19496 100644 --- a/src/USER-MISC/fix_gle.cpp +++ b/src/USER-MISC/fix_gle.cpp @@ -843,7 +843,7 @@ void FixGLE::unpack_restart(int nlocal, int nth) // skip to the nth set of extended variables // unpack the Nth first values this way because other fixes pack them - + int m = 0; for (int i = 0; i< nth; i++) m += static_cast (extra[nlocal][m]); m++; diff --git a/src/USER-MISC/fix_srp.cpp b/src/USER-MISC/fix_srp.cpp index 131fd7d27a..e72ea61b01 100644 --- a/src/USER-MISC/fix_srp.cpp +++ b/src/USER-MISC/fix_srp.cpp @@ -572,7 +572,7 @@ void FixSRP::unpack_restart(int nlocal, int nth) // skip to Nth set of extra values // unpack the Nth first values this way because other fixes pack them - + int m = 0; for (int i = 0; i < nth; i++){ m += static_cast (extra[nlocal][m]); diff --git a/src/USER-MISC/pair_meam_spline.cpp b/src/USER-MISC/pair_meam_spline.cpp index 5f35070130..5764b09d8a 100644 --- a/src/USER-MISC/pair_meam_spline.cpp +++ b/src/USER-MISC/pair_meam_spline.cpp @@ -448,7 +448,7 @@ void PairMEAMSpline::coeff(int narg, char **arg) if (map[j] == i) count++; if (count != 1) error->all(FLERR,"Pair style meam/spline requires one atom type per element"); - } + } } #define MAXLINE 1024 diff --git a/src/USER-OMP/ewald_omp.cpp b/src/USER-OMP/ewald_omp.cpp index e023daf2db..ec8d708da2 100644 --- a/src/USER-OMP/ewald_omp.cpp +++ b/src/USER-OMP/ewald_omp.cpp @@ -438,14 +438,14 @@ void EwaldOMP::eik_dot_r_triclinic() #pragma omp parallel LMP_DEFAULT_NONE #endif { - + int i,ifrom,ito,k,l,m,n,ic,tid; double cstr1,sstr1; double sqk,clpm,slpm; double unitk_lamda[3]; loop_setup_thr(ifrom,ito,tid,nlocal,nthreads); - + double max_kvecs[3]; max_kvecs[0] = kxmax; max_kvecs[1] = kymax; diff --git a/src/fix_neigh_history.cpp b/src/fix_neigh_history.cpp index 933bcdc265..e3d0d2830a 100644 --- a/src/fix_neigh_history.cpp +++ b/src/fix_neigh_history.cpp @@ -853,7 +853,7 @@ int FixNeighHistory::pack_restart(int i, double *buf) memcpy(&buf[m],&valuepartner[i][dnum*n],dnumbytes); m += dnum; } - // pack buf[0] this way because other fixes unpack it + // pack buf[0] this way because other fixes unpack it buf[0] = m; return m; }