jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7881 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2012-03-02 15:41:18 +00:00
parent 9d253478b4
commit 4deb839d4b
5 changed files with 60 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
src/USER-OMP/ewald_omp.cpp
View File

@ -59,7 +59,8 @@ void EwaldOMP::compute(int eflag, int vflag)
// set energy/virial flags
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = eflag_global = vflag_global = eflag_atom = vflag_atom = 0;
else evflag = evflag_atom = eflag_global = vflag_global =
eflag_atom = vflag_atom = 0;
// extend size of per-atom arrays if necessary
@ -99,7 +100,7 @@ void EwaldOMP::compute(int eflag, int vflag)
#endif
{
int i,k,ifrom,ito,tid;
int i,j,k,ifrom,ito,tid;
int kx,ky,kz;
double cypz,sypz,exprl,expim,partial;
@ -128,16 +129,13 @@ void EwaldOMP::compute(int eflag, int vflag)
ek[i][1] += partial*eg[k][1];
ek[i][2] += partial*eg[k][2];
#if 0
if (eflag_atom || vflag_atom) {
partial_peratom = exprl*sfacrl_all[k] + expim*sfacim_all[k];
if (evflag_atom) {
const double partial_peratom = exprl*sfacrl_all[k] + expim*sfacim_all[k];
if (eflag_atom) eatom[i] += q[i]*ug[k]*partial_peratom;
if (vflag_atom)
for (j = 0; j < 6; j++)
vatom[i][j] += ug[k]*vg[k][j]*partial_peratom;
}
#endif
}
}
@ -150,7 +148,7 @@ void EwaldOMP::compute(int eflag, int vflag)
f[i][2] += fac*ek[i][2];
}
// energy if requested
// global energy
if (eflag_global) {
#if defined(_OPENMP)
@ -161,7 +159,7 @@ void EwaldOMP::compute(int eflag, int vflag)
sfacim_all[k]*sfacim_all[k]);
}
// virial if requested
// global virial
if (vflag_global) {
#if defined(_OPENMP)
@ -178,6 +176,23 @@ void EwaldOMP::compute(int eflag, int vflag)
}
}
// per-atom energy/virial
// energy includes self-energy correction
if (evflag_atom) {
if (eflag_atom) {
for (i = ifrom; i < ito; i++) {
eatom[i] -= g_ewald*q[i]*q[i]/MY_PIS + MY_PI2*q[i]*qsum /
(g_ewald*g_ewald*volume);
eatom[i] *= qscale;
}
}
if (vflag_atom)
for (i = ifrom; i < ito; i++)
for (j = 0; j < 6; j++) vatom[i][j] *= q[i]*qscale;
}
reduce_thr(this, eflag,vflag,thr);
} // end of omp parallel region
@ -196,25 +211,6 @@ void EwaldOMP::compute(int eflag, int vflag)
virial[5] = v5 * qscale;
}
#if 0
// per-atom energy/virial
// energy includes self-energy correction
if (eflag_atom || vflag_atom) {
if (eflag_atom) {
for (i = 0; i < nlocal; i++) {
eatom[i] -= g_ewald*q[i]*q[i]/MY_PIS + MY_PI2*q[i]*qsum /
(g_ewald*g_ewald*volume);
eatom[i] *= qscale;
}
}
if (vflag_atom)
for (i = 0; i < nlocal; i++)
for (j = 0; j < 6; j++) vatom[i][j] *= q[i]*qscale;
}
#endif
if (slabflag) slabcorr();
}
@ -288,8 +284,8 @@ void EwaldOMP::eik_dot_r()
// 1 = (k,l,0), 2 = (k,-l,0)
for (k = 1; k <= kmax; k++) {
for (l = 1; l <= kmax; l++) {
for (k = 1; k <= kxmax; k++) {
for (l = 1; l <= kymax; l++) {
sqk = (k*unitk[0] * k*unitk[0]) + (l*unitk[1] * l*unitk[1]);
if (sqk <= gsqmx) {
cstr1 = 0.0;
@ -312,8 +308,8 @@ void EwaldOMP::eik_dot_r()
// 1 = (0,l,m), 2 = (0,l,-m)
for (l = 1; l <= kmax; l++) {
for (m = 1; m <= kmax; m++) {
for (l = 1; l <= kymax; l++) {
for (m = 1; m <= kzmax; m++) {
sqk = (l*unitk[1] * l*unitk[1]) + (m*unitk[2] * m*unitk[2]);
if (sqk <= gsqmx) {
cstr1 = 0.0;
@ -336,8 +332,8 @@ void EwaldOMP::eik_dot_r()
// 1 = (k,0,m), 2 = (k,0,-m)
for (k = 1; k <= kmax; k++) {
for (m = 1; m <= kmax; m++) {
for (k = 1; k <= kxmax; k++) {
for (m = 1; m <= kzmax; m++) {
sqk = (k*unitk[0] * k*unitk[0]) + (m*unitk[2] * m*unitk[2]);
if (sqk <= gsqmx) {
cstr1 = 0.0;
@ -360,9 +356,9 @@ void EwaldOMP::eik_dot_r()
// 1 = (k,l,m), 2 = (k,-l,m), 3 = (k,l,-m), 4 = (k,-l,-m)
for (k = 1; k <= kmax; k++) {
for (l = 1; l <= kmax; l++) {
for (m = 1; m <= kmax; m++) {
for (k = 1; k <= kxmax; k++) {
for (l = 1; l <= kymax; l++) {
for (m = 1; m <= kzmax; m++) {
sqk = (k*unitk[0] * k*unitk[0]) + (l*unitk[1] * l*unitk[1]) +
(m*unitk[2] * m*unitk[2]);
if (sqk <= gsqmx) {

4
src/USER-OMP/pppm_omp.cpp
View File

@ -344,7 +344,9 @@ void PPPMOMP::make_rho()
#if defined(_OPENMP)
#pragma omp parallel default(none)
#endif
{
#if defined(_OPENMP)
// each thread works on a fixed chunk of atoms.
const int tid = omp_get_thread_num();
const int inum = nlocal;
@ -406,8 +408,8 @@ void PPPMOMP::make_rho()
if (nthreads > 1) {
data_reduce_fft(&(density_brick[nzlo_out][nylo_out][nxlo_out]),ngrid,nthreads,1,tid);
}
}
#endif
}
}
/* ----------------------------------------------------------------------

1
src/USER-OMP/thr_data.cpp
View File

@ -170,7 +170,6 @@ void ThrData::virial_fdotr_compute(double **x, int nlocal, int nghost, int nfirs
virial_pair[4] += _f[i][2]*x[i][0];
virial_pair[5] += _f[i][2]*x[i][1];
}
nall = nlocal + nghost;
for (int i = nlocal; i < nall; i++) {
virial_pair[0] += _f[i][0]*x[i][0];

32
src/USER-OMP/thr_data.h
View File

@ -67,28 +67,28 @@ class ThrData {
private:
double eng_vdwl; // non-bonded non-coulomb energy
double eng_coul; // non-bonded coulomb energy
double virial_pair[6]; // virial contribution from non-bonded
double *eatom_pair;
double **vatom_pair;
double eng_bond; // bond energy
double virial_bond[6]; // virial contribution from bonds
double *eatom_bond;
double **vatom_bond;
double eng_angle; // angle energy
double virial_angle[6]; // virial contribution from angles
double *eatom_angle;
double **vatom_angle;
double eng_dihed; // dihedral energy
double virial_dihed[6]; // virial contribution from dihedrals
double *eatom_dihed;
double **vatom_dihed;
double eng_imprp; // improper energy
double virial_imprp[6]; // virial contribution from impropers
double *eatom_imprp;
double **vatom_imprp;
double eng_kspce; // kspace energy
double virial_pair[6]; // virial contribution from non-bonded
double virial_bond[6]; // virial contribution from bonds
double virial_angle[6]; // virial contribution from angles
double virial_dihed[6]; // virial contribution from dihedrals
double virial_imprp[6]; // virial contribution from impropers
double virial_kspce[6]; // virial contribution from kspace
double *eatom_pair;
double *eatom_bond;
double *eatom_angle;
double *eatom_dihed;
double *eatom_imprp;
double *eatom_kspce;
double **vatom_pair;
double **vatom_bond;
double **vatom_angle;
double **vatom_dihed;
double **vatom_imprp;
double **vatom_kspce;
// per thread segments of various force or similar arrays
@ -113,7 +113,7 @@ class ThrData {
const int _tid;
public:
// compute global per thread virial contribution from per-thread force
// compute global per thread virial contribution from global forces and positions
void virial_fdotr_compute(double **, int, int, int);
double memory_usage();

43
src/USER-OMP/thr_omp.cpp
View File

@ -134,20 +134,7 @@ void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom,
}
}
#if 0 /* not supported (yet) */
if (thr_style & THR_KSPACE) {
if (eflag & 2) {
thr->eatom_kspce = eatom + tid*nall;
if (nall > 0)
memset(&(thr->eatom_kspce[0]),0,nall*sizeof(double));
}
if (vflag & 4) {
thr->vatom_kspce = vatom + tid*nall;
if (nall > 0)
memset(&(thr->vatom_kspce[0][0]),0,nall*6*sizeof(double));
}
}
#endif
// nothing to do for THR_KSPACE
}
/* ----------------------------------------------------------------------
@ -255,7 +242,12 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag,
#endif
{
if (tid < nproxy) {
// nothing to do for kspace?
// nothing to collect for kspace proxy threads
// just reset pair accumulators to 0.0.
if (eflag & 1) {
thr->eng_vdwl = 0.0;
thr->eng_coul = 0.0;
}
if (vflag & 3)
for (int i=0; i < 6; ++i) {
thr->virial_pair[i] = 0.0;
@ -403,26 +395,7 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag,
case THR_KSPACE|THR_PROXY: // fallthrough
case THR_KSPACE:
// nothing to do (for now)
#if 0
if (evflag) {
KSpace *kspace = lmp->force->kspace;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
kspace->energy += thr->eng_kspce;
for (int i=0; i < 6; ++i)
kspace->virial[i] += thr->virial_kspce[i];
}
if (eflag & 2) {
data_reduce_thr(&(kspace->eatom[0]), nall, nthreads, 1, tid);
}
if (vflag & 4) {
data_reduce_thr(&(kspace->vatom[0][0]), nall, nthreads, 6, tid);
}
}
#endif
// nothing to do
break;
default: