From f350cc4ec77b5691387f53f460e013e0acb83672 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 20 Oct 2014 21:17:46 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12637 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GPU/pppm_gpu.cpp | 17 ++++++++++------- src/KSPACE/ewald.cpp | 15 +++++++++------ src/KSPACE/ewald_disp.cpp | 3 +++ src/KSPACE/msm.cpp | 15 +++++++++------ src/KSPACE/msm_cg.cpp | 13 +++++++++++-- src/KSPACE/pppm.cpp | 15 +++++++++------ src/KSPACE/pppm_cg.cpp | 15 +++++++++------ src/KSPACE/pppm_disp.cpp | 15 +++++++++------ src/KSPACE/pppm_stagger.cpp | 15 +++++++++------ src/USER-CUDA/pppm_cuda.cpp | 13 +++++++------ src/USER-CUDA/pppm_old.cpp | 15 +++++++++------ src/USER-OMP/ewald_omp.cpp | 10 ++++++++++ src/USER-OMP/msm_cg_omp.cpp | 13 +++++++++++-- src/force.cpp | 10 +++++----- src/kspace.cpp | 21 ++++++++++++--------- src/pair_hybrid.cpp | 2 +- 16 files changed, 133 insertions(+), 74 deletions(-) diff --git a/src/GPU/pppm_gpu.cpp b/src/GPU/pppm_gpu.cpp index f27665e227..c026dca445 100644 --- a/src/GPU/pppm_gpu.cpp +++ b/src/GPU/pppm_gpu.cpp @@ -284,21 +284,24 @@ void PPPMGPU::compute(int eflag, int vflag) if (evflag_atom) fieldforce_peratom(); + // update qsum and qsqsum, if needed + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + // sum energy across procs and add in volume-dependent term - // reset qsum and qsqsum if atom count has changed if (eflag_global) { double energy_all; MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - if (atom->natoms != natoms_original) { - qsum_qsq(0); - natoms_original = atom->natoms; - } - energy *= 0.5*volume; - energy -= g_ewald*qsqsum/1.772453851 + + energy -= g_ewald*qsqsum/MY_PIS + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); energy *= qscale; } diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index 054bf343e4..115dc95f94 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -431,19 +431,22 @@ void Ewald::compute(int eflag, int vflag) if (slabflag != 2) f[i][2] += qscale * q[i]*ek[i][2]; } + // update qsum and qsqsum, if needed + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + // sum global energy across Kspace vevs and add in volume-dependent term - // reset qsum and qsqsum if atom count has changed if (eflag_global) { for (k = 0; k < kcount; k++) energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] + sfacim_all[k]*sfacim_all[k]); - if (qsum_update_flag || (atom->natoms != natoms_original)) { - qsum_qsq(0); - natoms_original = atom->natoms; - } - energy -= g_ewald*qsqsum/MY_PIS + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); energy *= qscale; diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp index a8a5bc6d1f..f7c2275aa3 100644 --- a/src/KSPACE/ewald_disp.cpp +++ b/src/KSPACE/ewald_disp.cpp @@ -246,6 +246,9 @@ void EwaldDisp::setup() error->all(FLERR,"KSpace accuracy too low"); } + if (qsum_update_flag) + error->all(FLERR,"Kspace style ewald/disp does not support dynamic charges"); + bigint natoms = atom->natoms; double err; int kxmax = 1; diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index 7a33c2a765..2e752c7a2f 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -566,8 +566,16 @@ void MSM::compute(int eflag, int vflag) if (evflag_atom) fieldforce_peratom(); + // update qsum and qsqsum, if needed + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + // sum global energy across procs and add in self-energy term - // reset qsum and qsqsum if atom count has changed const double qscale = qqrd2e * scale; @@ -576,11 +584,6 @@ void MSM::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - if (qsum_update_flag || (atom->natoms != natoms_original)) { - qsum_qsq(0); - natoms_original = atom->natoms; - } - double e_self = qsqsum*gamma(0.0)/cutoff; energy -= e_self; energy *= 0.5*qscale; diff --git a/src/KSPACE/msm_cg.cpp b/src/KSPACE/msm_cg.cpp index 5e351cfcf2..7a743cc2dc 100644 --- a/src/KSPACE/msm_cg.cpp +++ b/src/KSPACE/msm_cg.cpp @@ -248,7 +248,16 @@ void MSMCG::compute(int eflag, int vflag) if (evflag_atom) fieldforce_peratom(); - // total long-range energy + // update qsum and qsqsum, if needed + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + + // sum global energy across procs and add in self-energy term const double qscale = force->qqrd2e * scale; @@ -257,7 +266,7 @@ void MSMCG::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - double e_self = qsqsum*gamma(0.0)/cutoff; // Self-energy term + double e_self = qsqsum*gamma(0.0)/cutoff; energy -= e_self; energy *= 0.5*qscale; } diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index f0ec0b9db5..d6305d3a52 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -667,8 +667,16 @@ void PPPM::compute(int eflag, int vflag) if (evflag_atom) fieldforce_peratom(); + // update qsum and qsqsum, if needed + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + // sum global energy across procs and add in volume-dependent term - // reset qsum and qsqsum if atom count has changed const double qscale = qqrd2e * scale; @@ -677,11 +685,6 @@ void PPPM::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - if (qsum_update_flag || (atom->natoms != natoms_original)) { - qsum_qsq(0); - natoms_original = atom->natoms; - } - energy *= 0.5*volume; energy -= g_ewald*qsqsum/MY_PIS + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp index 3ebcb5d3a4..e457685802 100644 --- a/src/KSPACE/pppm_cg.cpp +++ b/src/KSPACE/pppm_cg.cpp @@ -208,8 +208,16 @@ void PPPMCG::compute(int eflag, int vflag) if (evflag_atom) fieldforce_peratom(); + // update qsum and qsqsum, if needed + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + // sum global energy across procs and add in volume-dependent term - // reset qsum and qsqsum if atom count has changed const double qscale = qqrd2e * scale; @@ -218,11 +226,6 @@ void PPPMCG::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - if (qsum_update_flag || (atom->natoms != natoms_original)) { - qsum_qsq(0); - natoms_original = atom->natoms; - } - energy *= 0.5*volume; energy -= g_ewald*qsqsum/MY_PIS + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index 78611360bd..b79f682ef9 100755 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -1121,8 +1121,16 @@ void PPPMDisp::compute(int eflag, int vflag) if (evflag_atom) fieldforce_none_peratom(); } + // update qsum and qsqsum, if needed + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + // sum energy across procs and add in volume-dependent term - // reset qsum and qsqsum if atom count has changed const double qscale = force->qqrd2e * scale; if (eflag_global) { @@ -1135,11 +1143,6 @@ void PPPMDisp::compute(int eflag, int vflag) energy_1 *= 0.5*volume; energy_6 *= 0.5*volume; - if (qsum_update_flag || (atom->natoms != natoms_original)) { - qsum_qsq(0); - natoms_original = atom->natoms; - } - energy_1 -= g_ewald*qsqsum/MY_PIS + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); energy_6 += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij + diff --git a/src/KSPACE/pppm_stagger.cpp b/src/KSPACE/pppm_stagger.cpp index bc1a1ffe24..34a06c0cd2 100755 --- a/src/KSPACE/pppm_stagger.cpp +++ b/src/KSPACE/pppm_stagger.cpp @@ -202,8 +202,16 @@ void PPPMStagger::compute(int eflag, int vflag) stagger += 1.0/float(nstagger); } + // update qsum and qsqsum, if needed + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + // sum global energy across procs and add in volume-dependent term - // reset qsum and qsqsum if atom count has changed const double qscale = qqrd2e * scale; @@ -212,11 +220,6 @@ void PPPMStagger::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - if (qsum_update_flag || (atom->natoms != natoms_original)) { - qsum_qsq(0); - natoms_original = atom->natoms; - } - energy *= 0.5*volume/float(nstagger); energy -= g_ewald*qsqsum/MY_PIS + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); diff --git a/src/USER-CUDA/pppm_cuda.cpp b/src/USER-CUDA/pppm_cuda.cpp index 57797497f3..a8325e6ce0 100644 --- a/src/USER-CUDA/pppm_cuda.cpp +++ b/src/USER-CUDA/pppm_cuda.cpp @@ -690,17 +690,18 @@ void PPPMCuda::compute(int eflag, int vflag) } // extend size of per-atom arrays if necessary + // and force update of device data, if charges can change - if ((cu_atom->update_nmax)||(old_nmax==0)) { + if ((cu_atom->update_nmax)||(old_nmax==0)||qsum_update_flag) { memory->destroy(part2grid); nmax = atom->nmax; memory->create(part2grid,nmax,3,"pppm:part2grid"); - delete cu_part2grid; - delete [] adev_data_array; - adev_data_array=new dev_array[1]; - cu_part2grid = new cCudaData ((int*)part2grid,adev_data_array, nmax,3); + delete cu_part2grid; + delete [] adev_data_array; + adev_data_array=new dev_array[1]; + cu_part2grid = new cCudaData ((int*)part2grid,adev_data_array, nmax,3); - pppm_device_update(&cuda->shared_data,cu_part2grid->dev_data(),atom->nlocal,atom->nmax); + pppm_device_update(&cuda->shared_data,cu_part2grid->dev_data(),atom->nlocal,atom->nmax); old_nmax=nmax; } if(cu_atom->update_nlocal) {pppm_update_nlocal(cu_atom->nlocal);} diff --git a/src/USER-CUDA/pppm_old.cpp b/src/USER-CUDA/pppm_old.cpp index c2d1cad053..80ee9c97ae 100755 --- a/src/USER-CUDA/pppm_old.cpp +++ b/src/USER-CUDA/pppm_old.cpp @@ -724,8 +724,16 @@ void PPPMOld::compute(int eflag, int vflag) if (evflag_atom) fieldforce_peratom(); + // update qsum and qsqsum, if needed + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + // sum global energy across procs and add in volume-dependent term - // reset qsum and qsqsum if atom count has changed const double qscale = qqrd2e * scale; @@ -734,11 +742,6 @@ void PPPMOld::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - if (atom->natoms != natoms_original) { - qsum_qsq(0); - natoms_original = atom->natoms; - } - energy *= 0.5*volume; energy -= g_ewald*qsqsum/MY_PIS + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); diff --git a/src/USER-OMP/ewald_omp.cpp b/src/USER-OMP/ewald_omp.cpp index 07881bf895..9d06f321a9 100644 --- a/src/USER-OMP/ewald_omp.cpp +++ b/src/USER-OMP/ewald_omp.cpp @@ -86,6 +86,16 @@ void EwaldOMP::compute(int eflag, int vflag) MPI_Allreduce(sfacrl,sfacrl_all,kcount,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(sfacim,sfacim_all,kcount,MPI_DOUBLE,MPI_SUM,world); + // update qsum and qsqsum, if needed + // (n.b. needs to be done outside of the multi-threaded region) + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + // K-space portion of electric field // double loop over K-vectors and local atoms diff --git a/src/USER-OMP/msm_cg_omp.cpp b/src/USER-OMP/msm_cg_omp.cpp index cf4128eb73..72bdd2f958 100644 --- a/src/USER-OMP/msm_cg_omp.cpp +++ b/src/USER-OMP/msm_cg_omp.cpp @@ -250,7 +250,16 @@ void MSMCGOMP::compute(int eflag, int vflag) if (evflag_atom) fieldforce_peratom(); - // total long-range energy + // update qsum and qsqsum, if needed + + if (eflag_global || eflag_atom) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { + qsum_qsq(0); + natoms_original = atom->natoms; + } + } + + // sum global energy across procs and add in self-energy term const double qscale = force->qqrd2e * scale; @@ -259,7 +268,7 @@ void MSMCGOMP::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - double e_self = qsqsum*gamma(0.0)/cutoff; // Self-energy term + double e_self = qsqsum*gamma(0.0)/cutoff; energy -= e_self; energy *= 0.5*qscale; } diff --git a/src/force.cpp b/src/force.cpp index c7f74f9aac..dc0b1bad60 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -906,8 +906,7 @@ FILE *Force::open_potential(const char *name) fp = fopen(name,"r"); if (fp) { potential_date(fp,name); - fclose(fp); - fp = fopen(name,"r"); + rewind(fp); return fp; } @@ -934,9 +933,10 @@ FILE *Force::open_potential(const char *name) strcat(newpath,pot); fp = fopen(newpath,"r"); - potential_date(fp,name); - fclose(fp); - fp = fopen(newpath,"r"); + if (fp) { + potential_date(fp,name); + rewind(fp); + } delete [] newpath; return fp; diff --git a/src/kspace.cpp b/src/kspace.cpp index 7727c24c7a..2e7f9a8736 100644 --- a/src/kspace.cpp +++ b/src/kspace.cpp @@ -264,17 +264,20 @@ void KSpace::ev_setup(int eflag, int vflag) void KSpace::qsum_qsq(int flag) { - qsum = qsqsum = 0.0; - for (int i = 0; i < atom->nlocal; i++) { - qsum += atom->q[i]; - qsqsum += atom->q[i]*atom->q[i]; + const double * const q = atom->q; + const int nlocal = atom->nlocal; + double qsum_local(0.0), qsqsum_local(0.0); + +#if defined(_OPENMP) +#pragma omp parallel for default(none) reduction(+:qsum_local,qsqsum_local) +#endif + for (int i = 0; i < nlocal; i++) { + qsum_local += q[i]; + qsqsum_local += q[i]*q[i]; } - double tmp; - MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world); - qsum = tmp; - MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world); - qsqsum = tmp; + MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&qsqsum_local,&qsqsum,1,MPI_DOUBLE,MPI_SUM,world); if (qsqsum == 0.0) error->all(FLERR,"Cannot use kspace solver on system with no charge"); diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index 164dd52a21..3d2fa67d9e 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -745,7 +745,7 @@ void *PairHybrid::extract(const char *str, int &dim) { void *cutptr = NULL; void *ptr; - double cutvalue; + double cutvalue = 0.0; for (int m = 0; m < nstyles; m++) { ptr = styles[m]->extract(str,dim);