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@8153 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2012-05-23 14:43:28 +00:00
parent c1e5db9ff0
commit 365218d56f
2 changed files with 18 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
src/KSPACE/ewald.cpp
View File

@ -1023,7 +1023,7 @@ void Ewald::compute_group_group(int groupbit_A, int groupbit_B, int BA_flag)
int kx,ky,kz;
double cypz,sypz,exprl,expim;
// partial structure factors for groups A and B on each processor
for (k = 0; k < kcount; k++) {
@ -1071,46 +1071,17 @@ void Ewald::compute_group_group(int groupbit_A, int groupbit_B, int BA_flag)
const double qscale = force->qqrd2e * scale;
double partial_group;
// total group A <--> group B energy
// self and boundary correction terms are in compute_group_group.cpp
for (k = 0; k < kcount; k++) {
partial_group = sfacrl_A_all[k]*sfacrl_B_all[k] +
sfacim_A_all[k]*sfacim_B_all[k];
e2group += ug[k]*partial_group;
}
// total charge of groups A & B, needed for self-energy correction
double qsum_group = 0.0, qsqsum_group = 0.0;
for (int i = 0; i < atom->nlocal; i++) {
if ((mask[i] & groupbit_A) && (mask[i] & groupbit_B)) {
if (BA_flag) continue;
qsqsum_group += q[i]*q[i];
}
//if ((mask[i] & groupbit_A) || (mask[i] & groupbit_B))
// qsum_group += q[i];
}
double tmp;
MPI_Allreduce(&qsum_group,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsum_group = tmp;
//MPI_Allreduce(&qsqsum_group,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
//qsqsum_group = tmp;
// self-energy correction
e2group -= g_ewald*qsqsum_group/MY_PIS;
// should we include the boundary condition or not?
// e2group += MY_PI2*qsum_group*qsum_group / (g_ewald*g_ewald*volume);
e2group *= qscale;
// total group A <--> group B force
for (k = 0; k < kcount; k++) {

64
src/KSPACE/pppm.cpp
View File

@ -1020,7 +1020,7 @@ void PPPM::set_grid()
nx_pppm = static_cast<int> (xprd/h_x) + 1;
ny_pppm = static_cast<int> (yprd/h_y) + 1;
nz_pppm = static_cast<int> (zprd_slab/h_z) + 1;
err = rms(h_x,xprd,natoms,q2,acons);
while (err > accuracy) {
err = rms(h_x,xprd,natoms,q2,acons);
@ -2577,14 +2577,7 @@ void PPPM::compute_group_group(int groupbit_A, int groupbit_B, int BA_flag)
int nlocal = atom->nlocal;
int *mask = atom->mask;
// convert atoms from box to lamda coords
if (triclinic == 0) boxlo = domain->boxlo;
else {
boxlo = domain->boxlo_lamda;
domain->x2lamda(atom->nlocal);
}
// map my particle charge onto my local 3d density grid
make_rho_groups(groupbit_A,groupbit_B,BA_flag);
@ -2627,52 +2620,25 @@ void PPPM::compute_group_group(int groupbit_A, int groupbit_B, int BA_flag)
const double qscale = force->qqrd2e * scale;
// total group A <--> group B energy
// self and boundary correction terms are in compute_group_group.cpp
double e2group_all;
MPI_Allreduce(&e2group,&e2group_all,1,MPI_DOUBLE,MPI_SUM,world);
e2group = e2group_all;
e2group *= 0.5*volume;
// total charge of groups A & B, needed for self-energy correction
double qsum_group = 0.0, qsqsum_group = 0.0;
for (int i = 0; i < atom->nlocal; i++) {
if ((mask[i] & groupbit_A) && (mask[i] & groupbit_B)) {
if (BA_flag) continue;
qsqsum_group += q[i]*q[i];
}
}
double tmp;
MPI_Allreduce(&qsum_group,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsum_group = tmp;
// self-energy correction
e2group -= g_ewald*qsqsum_group/MY_PIS;
// should we include this boundary condition or not?
// e2group += MY_PI2*qsum_group*qsum_group / (g_ewald*g_ewald*volume);
e2group *= qscale;
e2group *= qscale*0.5*volume;
// total group A <--> group B force
double f2group_all[3];
MPI_Allreduce(f2group,f2group_all,3,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 3; i++) f2group[i] = qscale*volume*f2group_all[i];
// convert atoms back from lamda to box coords
if (triclinic) domain->lamda2x(atom->nlocal);
for (i = 0; i < 3; i++) f2group[i] = qscale*volume*f2group_all[i];
}
/* ----------------------------------------------------------------------
allocate group-group memory that depends on # of K-vectors and order
------------------------------------------------------------------------- */
allocate group-group memory that depends on # of K-vectors and order
------------------------------------------------------------------------- */
void PPPM::allocate_groups()
{
@ -2685,8 +2651,8 @@ void PPPM::allocate_groups()
}
/* ----------------------------------------------------------------------
deallocate group-group memory that depends on # of K-vectors and order
------------------------------------------------------------------------- */
deallocate group-group memory that depends on # of K-vectors and order
------------------------------------------------------------------------- */
void PPPM::deallocate_groups()
{
@ -2697,11 +2663,11 @@ void PPPM::deallocate_groups()
}
/* ----------------------------------------------------------------------
create discretized "density" on section of global grid due to my particles
density(x,y,z) = charge "density" at grid points of my 3d brick
(nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts)
in global grid for group-group interactions
------------------------------------------------------------------------- */
create discretized "density" on section of global grid due to my particles
density(x,y,z) = charge "density" at grid points of my 3d brick
(nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts)
in global grid for group-group interactions
------------------------------------------------------------------------- */
void PPPM::make_rho_groups(int groupbit_A, int groupbit_B, int BA_flag)
{
@ -2770,7 +2736,7 @@ void PPPM::make_rho_groups(int groupbit_A, int groupbit_B, int BA_flag)
/* ----------------------------------------------------------------------
FFT-based Poisson solver for group-group interactions
------------------------------------------------------------------------- */
------------------------------------------------------------------------- */
void PPPM::poisson_groups(int BA_flag)
{