diff --git a/src/GRANULAR/Install.csh b/src/GRANULAR/Install.csh index 4eb06edcea..8a5d4efefa 100644 --- a/src/GRANULAR/Install.csh +++ b/src/GRANULAR/Install.csh @@ -6,7 +6,6 @@ if ($1 == 1) then cp atom_vec_granular.cpp .. cp fix_freeze.cpp .. - cp fix_gran_diag.cpp .. cp fix_nve_gran.cpp .. cp fix_pour.cpp .. cp fix_wall_gran.cpp .. @@ -16,7 +15,6 @@ if ($1 == 1) then cp atom_vec_granular.h .. cp fix_freeze.h .. - cp fix_gran_diag.h .. cp fix_nve_gran.h .. cp fix_pour.h .. cp fix_wall_gran.h .. @@ -31,7 +29,6 @@ else if ($1 == 0) then rm ../atom_vec_granular.cpp rm ../fix_freeze.cpp - rm ../fix_gran_diag.cpp rm ../fix_nve_gran.cpp rm ../fix_pour.cpp rm ../fix_wall_gran.cpp @@ -41,7 +38,6 @@ else if ($1 == 0) then rm ../atom_vec_granular.h rm ../fix_freeze.h - rm ../fix_gran_diag.h rm ../fix_nve_gran.h rm ../fix_pour.h rm ../fix_wall_gran.h diff --git a/src/GRANULAR/fix_gran_diag.cpp b/src/GRANULAR/fix_gran_diag.cpp deleted file mode 100644 index d0c22bb200..0000000000 --- a/src/GRANULAR/fix_gran_diag.cpp +++ /dev/null @@ -1,1009 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Leo Silbert (SNL), Gary Grest (SNL) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdlib.h" -#include "string.h" -#include "fix_gran_diag.h" -#include "update.h" -#include "force.h" -#include "pair.h" -#include "atom.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "modify.h" -#include "comm.h" -#include "domain.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define MIN(a,b) ((a) < (b) ? (a) : (b)) -#define MAX(a,b) ((a) > (b) ? (a) : (b)) - -enum{NO_HISTORY,HISTORY,HERTZIAN}; - -/* ---------------------------------------------------------------------- */ - -FixGranDiag::FixGranDiag(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg != 6) error->all("Illegal fix gran/diag command"); - nevery = atoi(arg[3]); - if (nevery <= 0) error->all("Illegal fix gran/diag command"); - first = 1; - - if (!atom->radius_flag || !atom->rmass_flag || !atom->omega_flag) - error->all("Fix gran/diag requires atom attributes radius, rmass, omega"); - - MPI_Comm_rank(world,&me); - if (me == 0) { - char *file = new char[128]; - - sprintf(file,"%s.den",arg[4]); - fpden = fopen(file,"w"); - if (fpden == NULL) { - char str[128]; - sprintf(str,"Cannot open fix gran/diag file %s",file); - error->one(str); - } - - sprintf(file,"%s.vel",arg[4]); - fpvel = fopen(file,"w"); - if (fpvel == NULL) { - char str[128]; - sprintf(str,"Cannot open fix gran/diag file %s",file); - error->one(str); - } - - sprintf(file,"%s.str",arg[4]); - fpstr = fopen(file,"w"); - if (fpstr == NULL) { - char str[128]; - sprintf(str,"Cannot open fix gran/diag file %s",file); - error->one(str); - } - - delete [] file; - } - - step = atof(arg[5]); - stepinv = 1.0/step; - PI = 4.0*atan(1.0); - - maxlayers = 0; -} - -/* ---------------------------------------------------------------------- */ - -FixGranDiag::~FixGranDiag() -{ - deallocate(); - - if (me == 0) { - fclose(fpden); - fclose(fpvel); - fclose(fpstr); - } -} - -/* ---------------------------------------------------------------------- */ - -int FixGranDiag::setmask() -{ - int mask = 0; - mask |= END_OF_STEP; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixGranDiag::init() -{ - // set local values from Pair values - - if (force->pair == NULL) - error->all("Fix gran/diag is incompatible with Pair style"); - double *p_xkk = (double *) force->pair->extract("xkk"); - double *p_gamman = (double *) force->pair->extract("gamman"); - double *p_xmu = (double *) force->pair->extract("xmu"); - int *p_dampflag = (int *) force->pair->extract("dampflag"); - if (!p_xkk || !p_gamman || !p_xmu || !p_dampflag) - error->all("Fix gran/diag is incompatible with Pair style"); - xkk = *p_xkk; - double gamman = *p_gamman; - xmu = *p_xmu; - int dampflag = *p_dampflag; - - // same initialization as in pair_gran_history::init_style() - - xkkt = xkk * 2.0/7.0; - dt = update->dt; - double gammas = 0.5*gamman; - if (dampflag == 0) gammas = 0.0; - gamman_dl = gamman/dt; - gammas_dl = gammas/dt; - - // set pairstyle from granular pair style - - if (force->pair_match("gran/no_history")) pairstyle = NO_HISTORY; - else if (force->pair_match("gran/history")) pairstyle = HISTORY; - else if (force->pair_match("gran/hertzian")) pairstyle = HERTZIAN; - - // check for freeze Fix and set freeze_group_bit - - int i; - for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"freeze") == 0) break; - if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit; - else freeze_group_bit = 0; -} - -/* ---------------------------------------------------------------------- */ - -void FixGranDiag::setup() -{ - // get neighbor list from Pair class - // can't do this until setup since lists are setup by neighbor::init() - - list = pair->list; - - if (first) end_of_step(); - first = 0; -} - -/* ---------------------------------------------------------------------- */ - -void FixGranDiag::end_of_step() -{ - int i,m; - - // setup for 2d vs 3d - - int dim; - if (domain->dimension == 3) dim = 2; - else dim = 1; - - // set bottom of box for binning purposes - - boxlo = domain->boxlo[dim]; - - // update ghost atom info - // else ghost x/v is out-of-date at end of timestep - - comm->communicate(); - - // insure accumulator arrays are correct length - // add 10 for buffer - - nlayers = static_cast (domain->prd[dim] / step) + 10; - if (nlayers > maxlayers) { - deallocate(); - maxlayers = nlayers; - allocate(); - } - - // zero all accumulators - - for (i = 0; i < nlayers; i++) { - numdens[i] = 0; - dendens[i] = 0.0; - velx[i] = vely[i] = velz[i] = 0.0; - velxx[i] = velyy[i] = velzz[i] = velxy[i] = velxz[i] = velyz[i] = 0.0; - sigxx[i] = sigyy[i] = sigzz[i] = sigxy[i] = sigxz[i] = sigyz[i] = 0.0; - } - - // density/velocity accumulation by atom - // assign to layer based on atom distance from bottom of domain - - int overflow = 0; - - double **x = atom->x; - double **v = atom->v; - double *rmass = atom->rmass; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - m = static_cast ((x[i][dim]-boxlo) * stepinv); - if (m >= 0 && m < nlayers) { - numdens[m]++; - dendens[m] += rmass[i]; - velx[m] += v[i][0]; - vely[m] += v[i][1]; - velz[m] += v[i][2]; - velxx[m] += v[i][0]*v[i][0]; - velyy[m] += v[i][1]*v[i][1]; - velzz[m] += v[i][2]*v[i][2]; - velxy[m] += v[i][0]*v[i][1]; - velxz[m] += v[i][0]*v[i][2]; - velyz[m] += v[i][1]*v[i][2]; - } else overflow++; - } - - // m = largest layer # with any counts - // nmax = # of layers up to m - - for (m = nlayers-1; m >= 0; m--) if (numdens[m]) break; - int nmax = m + 1; - - int tmp = nmax; - MPI_Allreduce(&tmp,&nmax,1,MPI_INT,MPI_MAX,world); - - // overflow = total # of atoms out-of-bounds of layer arrays - - tmp = overflow; - MPI_Allreduce(&tmp,&overflow,1,MPI_INT,MPI_SUM,world); - - // sum contributions across procs - - int *isum = new int[nmax]; - double *dsum = new double[nmax]; - - MPI_Allreduce(numdens,isum,nmax,MPI_INT,MPI_SUM,world); - for (i = 0; i < nmax; i++) numdens[i] = isum[i]; - - MPI_Allreduce(dendens,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) dendens[i] = dsum[i]; - - MPI_Allreduce(velx,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) velx[i] = dsum[i]; - MPI_Allreduce(vely,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) vely[i] = dsum[i]; - MPI_Allreduce(velz,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) velz[i] = dsum[i]; - - MPI_Allreduce(velxx,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) velxx[i] = dsum[i]; - MPI_Allreduce(velyy,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) velyy[i] = dsum[i]; - MPI_Allreduce(velzz,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) velzz[i] = dsum[i]; - MPI_Allreduce(velxy,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) velxy[i] = dsum[i]; - MPI_Allreduce(velxz,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) velxz[i] = dsum[i]; - MPI_Allreduce(velyz,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) velyz[i] = dsum[i]; - - // compute contribution to stress by every atom pair - - if (pairstyle == NO_HISTORY) stress_no_history(); - else if (pairstyle == HISTORY) stress_history(); - else if (pairstyle == HERTZIAN) stress_hertzian(); - - // sum contributions across procs - - MPI_Allreduce(sigxx,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) sigxx[i] = dsum[i]; - MPI_Allreduce(sigyy,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) sigyy[i] = dsum[i]; - MPI_Allreduce(sigzz,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) sigzz[i] = dsum[i]; - - MPI_Allreduce(sigxy,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) sigxy[i] = dsum[i]; - MPI_Allreduce(sigxz,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) sigxz[i] = dsum[i]; - MPI_Allreduce(sigyz,dsum,nmax,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmax; i++) sigyz[i] = dsum[i]; - - delete [] isum; - delete [] dsum; - - // density/velocity/stress by layer - - double velxxd1,velyyd1,velzzd1,velxyd1,velxzd1,velyzd1; - double vollayer = domain->xprd; - if (domain->dimension == 3) vollayer *= domain->yprd; - - for (m = 0; m < nmax; m++) { - if (numdens[m] == 0) numdens[m] = 1; - dendens[m] = stepinv*dendens[m]/vollayer; - - velx11[m] = velx[m]/numdens[m]; - vely11[m] = vely[m]/numdens[m]; - velz11[m] = velz[m]/numdens[m]; - velxx11[m] = velxx[m]/numdens[m]; - velyy11[m] = velyy[m]/numdens[m]; - velzz11[m] = velzz[m]/numdens[m]; - velxy11[m] = velxy[m]/numdens[m]; - velxz11[m] = velxz[m]/numdens[m]; - velyz11[m] = velyz[m]/numdens[m]; - - velxxd1 = velxx11[m] - velx11[m]*velx11[m]; - velyyd1 = velyy11[m] - vely11[m]*vely11[m]; - velzzd1 = velzz11[m] - velz11[m]*velz11[m]; - velxyd1 = velxy11[m] - velx11[m]*vely11[m]; - velxzd1 = velxz11[m] - velx11[m]*velz11[m]; - velyzd1 = velyz11[m] - vely11[m]*velz11[m]; - - velfxx[m] = velxxd1 * dendens[m]; - velfyy[m] = velyyd1 * dendens[m]; - velfzz[m] = velzzd1 * dendens[m]; - velfxy[m] = velxyd1 * dendens[m]; - velfxz[m] = velxzd1 * dendens[m]; - velfyz[m] = velyzd1 * dendens[m]; - - sigx2[m] = sigxx[m]/(2.0*vollayer*step) + velxxd1*dendens[m]; - sigy2[m] = sigyy[m]/(2.0*vollayer*step) + velyyd1*dendens[m]; - sigz2[m] = sigzz[m]/(2.0*vollayer*step) + velzzd1*dendens[m]; - sigxy2[m] = sigxy[m]/(2.0*vollayer*step) + velxyd1*dendens[m]; - sigxz2[m] = sigxz[m]/(2.0*vollayer*step) + velxzd1*dendens[m]; - sigyz2[m] = sigyz[m]/(2.0*vollayer*step) + velyzd1*dendens[m]; - } - - // write out density profile - - if (me == 0) { - fprintf(fpden,"ITEM: TIMESTEP\n"); - fprintf(fpden,"%d\n",update->ntimestep); - fprintf(fpden,"ITEM: NUMBER OF LAYERS / OVERFLOWS\n"); - fprintf(fpden,"%d %d\n",nmax,overflow); - fprintf(fpden,"ITEM: DENSITY BY LAYER\n"); - for (m = 0; m < nmax; m++) - fprintf(fpden,"%d %g %g\n",m+1,(m+1)*step+boxlo,dendens[m]*PI/6.0); - } - - // write out velocity profile - - if (me == 0) { - fprintf(fpvel,"ITEM: TIMESTEP\n"); - fprintf(fpvel,"%d\n",update->ntimestep); - fprintf(fpvel,"ITEM: NUMBER OF LAYERS / OVERFLOWS\n"); - fprintf(fpvel,"%d %d\n",nmax,overflow); - fprintf(fpvel,"ITEM: VELOCITY BY LAYER\n"); - for (m = 0; m < nmax; m++) - fprintf(fpvel,"%d %g %g %g %g %g %g %g %g %g %g\n", - m+1,(m+1)*step+boxlo, - velx11[m],vely11[m],velz11[m], - velxx11[m],velyy11[m],velzz11[m], - velxy11[m],velxz11[m],velyz11[m]); - } - - // write out stress profile - - if (me == 0) { - fprintf(fpstr,"ITEM: TIMESTEP\n"); - fprintf(fpstr,"%d\n",update->ntimestep); - fprintf(fpstr,"ITEM: NUMBER OF LAYERS / OVERFLOWS\n"); - fprintf(fpstr,"%d %d\n",nmax,overflow); - fprintf(fpstr,"ITEM: STRESS BY LAYER\n"); - for (m = 0; m < nmax; m++) - fprintf(fpstr,"%d %g %g %g %g %g %g %g %g %g %g %g %g %g\n", - m+1,(m+1)*step+boxlo, - sigx2[m],sigy2[m],sigz2[m], - sigxy2[m],sigxz2[m],sigyz2[m], - velfxx[m],velfyy[m],velfzz[m], - velfxy[m],velfxz[m],velfyz[m]); - } -} - -/* ---------------------------------------------------------------------- */ - -void FixGranDiag::allocate() -{ - numdens = new int[maxlayers]; - dendens = new double[maxlayers]; - - velx = new double[maxlayers]; - vely = new double[maxlayers]; - velz = new double[maxlayers]; - velxx = new double[maxlayers]; - velyy = new double[maxlayers]; - velzz = new double[maxlayers]; - velxy = new double[maxlayers]; - velxz = new double[maxlayers]; - velyz = new double[maxlayers]; - - velx11 = new double[maxlayers]; - vely11 = new double[maxlayers]; - velz11 = new double[maxlayers]; - velxx11 = new double[maxlayers]; - velyy11 = new double[maxlayers]; - velzz11 = new double[maxlayers]; - velxy11 = new double[maxlayers]; - velxz11 = new double[maxlayers]; - velyz11 = new double[maxlayers]; - - sigxx = new double[maxlayers]; - sigyy = new double[maxlayers]; - sigzz = new double[maxlayers]; - sigxy = new double[maxlayers]; - sigxz = new double[maxlayers]; - sigyz = new double[maxlayers]; - - sigx2 = new double[maxlayers]; - sigy2 = new double[maxlayers]; - sigz2 = new double[maxlayers]; - sigxy2 = new double[maxlayers]; - sigxz2 = new double[maxlayers]; - sigyz2 = new double[maxlayers]; - - velfxx = new double[maxlayers]; - velfyy = new double[maxlayers]; - velfzz = new double[maxlayers]; - velfxy = new double[maxlayers]; - velfxz = new double[maxlayers]; - velfyz = new double[maxlayers]; -} - -/* ---------------------------------------------------------------------- */ - -void FixGranDiag::deallocate() -{ - if (maxlayers == 0) return; - - delete [] numdens; delete [] dendens; - - delete [] velx; delete [] vely; delete [] velz; - delete [] velxx; delete [] velyy; delete [] velzz; - delete [] velxy; delete [] velxz; delete [] velyz; - - delete [] velx11; delete [] vely11; delete [] velz11; - delete [] velxx11; delete [] velyy11; delete [] velzz11; - delete [] velxy11; delete [] velxz11; delete [] velyz11; - - delete [] sigxx; delete [] sigyy; delete [] sigzz; - delete [] sigxy; delete [] sigxz; delete [] sigyz; - - delete [] sigx2; delete [] sigy2; delete [] sigz2; - delete [] sigxy2; delete [] sigxz2; delete [] sigyz2; - - delete [] velfxx; delete [] velfyy; delete [] velfzz; - delete [] velfxy; delete [] velfxz; delete [] velfyz; -} - -/* ---------------------------------------------------------------------- */ - -void FixGranDiag::stress_no_history() -{ - int i,j,ii,jj,m,inum,jnum; - double xtmp,ytmp,ztmp,delx,dely,delz; - double radi,radj,radsum,rsq,r; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - double xmeff,damp,ccel,ccelx,ccely,ccelz; - double fn,fs,ft,fs1,fs2,fs3; - int *ilist,*jlist,*numneigh,**firstneigh; - - double **x = atom->x; - double **v = atom->v; - double **omega = atom->omega; - double *radius = atom->radius; - double *rmass = atom->rmass; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // loop over all neighbors of my atoms - // store stress for both atoms i and j - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - radi = radius[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - - // skip if neither atom is in fix group - - if (!(mask[i] & groupbit) && !(mask[j] & groupbit)) continue; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radj = radius[j]; - radsum = radi + radj; - - if (rsq < radsum*radsum) { - - r = sqrt(rsq); - - // relative translational velocity - - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - vr1 *= dt; - vr2 *= dt; - vr3 *= dt; - - // normal component - - vnnr = vr1*delx + vr2*dely + vr3*delz; - vn1 = delx*vnnr / rsq; - vn2 = dely*vnnr / rsq; - vn3 = delz*vnnr / rsq; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - wr1 = radi*omega[i][0] + radj*omega[j][0]; - wr2 = radi*omega[i][1] + radj*omega[j][1]; - wr3 = radi*omega[i][2] + radj*omega[j][2]; - - wr1 *= dt/r; - wr2 *= dt/r; - wr3 *= dt/r; - - // normal damping term - // this definition of DAMP includes the extra 1/r term - - xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); - if (mask[i] & freeze_group_bit) xmeff = rmass[j]; - if (mask[j] & freeze_group_bit) xmeff = rmass[i]; - damp = xmeff*gamman_dl*vnnr/rsq; - ccel = xkk*(radsum-r)/r - damp; - - // relative velocities - - vtr1 = vt1 - (delz*wr2-dely*wr3); - vtr2 = vt2 - (delx*wr3-delz*wr1); - vtr3 = vt3 - (dely*wr1-delx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // force normalization - - fn = xmu * fabs(ccel*r); - fs = xmeff*gammas_dl*vrel; - if (vrel != 0.0) ft = MIN(fn,fs) / vrel; - else ft = 0.0; - - // shear friction forces - - fs1 = -ft*vtr1; - fs2 = -ft*vtr2; - fs3 = -ft*vtr3; - - // forces - - ccelx = delx*ccel + fs1; - ccely = dely*ccel + fs2; - ccelz = delz*ccel + fs3; - - // stress contribution of atom pair to z-layers - // atom i always contributes - // atom j contributes if newton_pair is on or if owned by this proc - - m = static_cast ((x[i][dim]-boxlo) * stepinv); - if (m >= 0 && m < nlayers) { - sigxx[m] += delx*ccelx; - sigyy[m] += dely*ccely; - sigzz[m] += delz*ccelz; - sigxy[m] += delx*ccely; - sigxz[m] += delx*ccelz; - sigyz[m] += dely*ccelz; - } - - if (newton_pair || j < nlocal) { - m = static_cast ((x[j][dim]-boxlo) * stepinv); - if (m >= 0 && m < nlayers) { - sigxx[m] += delx*ccelx; - sigyy[m] += dely*ccely; - sigzz[m] += delz*ccelz; - sigxy[m] += delx*ccely; - sigxz[m] += delx*ccelz; - sigyz[m] += dely*ccelz; - } - } - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixGranDiag::stress_history() -{ - int i,j,ii,jj,m,inum,jnum; - double xtmp,ytmp,ztmp,delx,dely,delz; - double radi,radj,radsum,rsq,r; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel,shrx,shry,shrz; - double xmeff,damp,ccel,ccelx,ccely,ccelz; - double fn,fs,fs1,fs2,fs3; - double shrmag,rsht; - int *ilist,*jlist,*numneigh,**firstneigh; - double *shear,*allshear,**firstshear; - - double **x = atom->x; - double **v = atom->v; - double **omega = atom->omega; - double *radius = atom->radius; - double *rmass = atom->rmass; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->listgranhistory->firstneigh; - firstshear = list->listgranhistory->firstdouble; - - // loop over all neighbors of my atoms - // store stress on both atoms i and j - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - radi = radius[i]; - allshear = firstshear[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - - // skip if neither atom is in fix group - - if (!(mask[i] & groupbit) && !(mask[j] & groupbit)) continue; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radj = radius[j]; - radsum = radi + radj; - - if (rsq < radsum*radsum) { - - r = sqrt(rsq); - - // relative translational velocity - - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - vr1 *= dt; - vr2 *= dt; - vr3 *= dt; - - // normal component - - vnnr = vr1*delx + vr2*dely + vr3*delz; - vn1 = delx*vnnr / rsq; - vn2 = dely*vnnr / rsq; - vn3 = delz*vnnr / rsq; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - wr1 = radi*omega[i][0] + radj*omega[j][0]; - wr2 = radi*omega[i][1] + radj*omega[j][1]; - wr3 = radi*omega[i][2] + radj*omega[j][2]; - - wr1 *= dt/r; - wr2 *= dt/r; - wr3 *= dt/r; - - // normal damping term - // this definition of DAMP includes the extra 1/r term - - xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); - if (mask[i] & freeze_group_bit) xmeff = rmass[j]; - if (mask[j] & freeze_group_bit) xmeff = rmass[i]; - damp = xmeff*gamman_dl*vnnr/rsq; - ccel = xkk*(radsum-r)/r - damp; - - // relative velocities - - vtr1 = vt1 - (delz*wr2-dely*wr3); - vtr2 = vt2 - (delx*wr3-delz*wr1); - vtr3 = vt3 - (dely*wr1-delx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - // shrmag = magnitude of shear - // do not update shear history since not timestepping - - shear = &allshear[3*jj]; - shrx = shear[0] + vtr1; - shry = shear[1] + vtr2; - shrz = shear[2] + vtr3; - shrmag = sqrt(shrx*shrx + shry*shry + shrz*shrz); - - // rotate shear displacements correctly - - rsht = shrx*delx + shry*dely + shrz*delz; - rsht /= rsq; - shrx -= rsht*delx; - shry -= rsht*dely; - shrz -= rsht*delz; - - // tangential forces - - fs1 = - (xkkt*shrx + xmeff*gammas_dl*vtr1); - fs2 = - (xkkt*shry + xmeff*gammas_dl*vtr2); - fs3 = - (xkkt*shrz + xmeff*gammas_dl*vtr3); - - // force normalization - // rescale frictional displacements and forces if needed - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - fn = xmu * fabs(ccel*r); - - if (fs > fn) { - if (shrmag != 0.0) { - fs1 *= fn/fs; - fs2 *= fn/fs; - fs3 *= fn/fs; - } else { - fs1 = 0.0; - fs2 = 0.0; - fs3 = 0.0; - } - } - - // forces - - ccelx = delx*ccel + fs1; - ccely = dely*ccel + fs2; - ccelz = delz*ccel + fs3; - - // stress contribution of atom pair to z-layers - // atom i always contributes - // atom j contributes if newton_pair is on or if owned by this proc - - m = static_cast ((x[i][dim]-boxlo) * stepinv); - if (m >= 0 && m < nlayers) { - sigxx[m] += delx*ccelx; - sigyy[m] += dely*ccely; - sigzz[m] += delz*ccelz; - sigxy[m] += delx*ccely; - sigxz[m] += delx*ccelz; - sigyz[m] += dely*ccelz; - } - - if (newton_pair || j < nlocal) { - m = static_cast ((x[j][dim]-boxlo) * stepinv); - if (m >= 0 && m < nlayers) { - sigxx[m] += delx*ccelx; - sigyy[m] += dely*ccely; - sigzz[m] += delz*ccelz; - sigxy[m] += delx*ccely; - sigxz[m] += delx*ccelz; - sigyz[m] += dely*ccelz; - } - } - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixGranDiag::stress_hertzian() -{ - int i,j,ii,jj,m,inum,jnum; - double xtmp,ytmp,ztmp,delx,dely,delz; - double radi,radj,radsum,rsq,r; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel,shrx,shry,shrz; - double xmeff,damp,ccel,ccelx,ccely,ccelz; - double fn,fs,fs1,fs2,fs3; - double shrmag,rsht,rhertz; - int *ilist,*jlist,*numneigh,**firstneigh; - double *shear,*allshear,**firstshear; - - double **x = atom->x; - double **v = atom->v; - double **omega = atom->omega; - double *radius = atom->radius; - double *rmass = atom->rmass; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - firstshear = listgran->firstdouble; - - // loop over all neighbors of my atoms - // store stress on both atoms i and j - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - radi = radius[i]; - allshear = firstshear[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - - // skip if neither atom is in fix group - - if (!(mask[i] & groupbit) && !(mask[j] & groupbit)) continue; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radj = radius[j]; - radsum = radi + radj; - - if (rsq < radsum*radsum) { - - r = sqrt(rsq); - - // relative translational velocity - - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - vr1 *= dt; - vr2 *= dt; - vr3 *= dt; - - // normal component - - vnnr = vr1*delx + vr2*dely + vr3*delz; - vn1 = delx*vnnr / rsq; - vn2 = dely*vnnr / rsq; - vn3 = delz*vnnr / rsq; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - wr1 = radi*omega[i][0] + radj*omega[j][0]; - wr2 = radi*omega[i][1] + radj*omega[j][1]; - wr3 = radi*omega[i][2] + radj*omega[j][2]; - - wr1 *= dt/r; - wr2 *= dt/r; - wr3 *= dt/r; - - // normal damping term - // this definition of DAMP includes the extra 1/r term - - xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); - if (mask[i] & freeze_group_bit) xmeff = rmass[j]; - if (mask[j] & freeze_group_bit) xmeff = rmass[i]; - damp = xmeff*gamman_dl*vnnr/rsq; - ccel = xkk*(radsum-r)/r - damp; - rhertz = sqrt(radsum - r); - ccel = rhertz * ccel; - - // relative velocities - - vtr1 = vt1 - (delz*wr2-dely*wr3); - vtr2 = vt2 - (delx*wr3-delz*wr1); - vtr3 = vt3 - (dely*wr1-delx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - // shrmag = magnitude of shear - // do not update shear history since not timestepping - - shear = &allshear[3*jj]; - shrx = shear[0] + vtr1; - shry = shear[1] + vtr2; - shrz = shear[2] + vtr3; - shrmag = sqrt(shrx*shrx + shry*shry + shrz*shrz); - - // rotate shear displacements correctly - - rsht = shrx*delx + shry*dely + shrz*delz; - rsht /= rsq; - shrx -= rsht*delx; - shry -= rsht*dely; - shrz -= rsht*delz; - - // tangential forces - - fs1 = -rhertz * (xkkt*shrx + xmeff*gammas_dl*vtr1); - fs2 = -rhertz * (xkkt*shry + xmeff*gammas_dl*vtr2); - fs3 = -rhertz * (xkkt*shrz + xmeff*gammas_dl*vtr3); - - // force normalization - // rescale frictional displacements and forces if needed - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - fn = xmu * fabs(ccel*r); - - if (fs > fn) { - if (shrmag != 0.0) { - fs1 *= fn/fs; - fs2 *= fn/fs; - fs3 *= fn/fs; - } else { - fs1 = 0.0; - fs2 = 0.0; - fs3 = 0.0; - } - } - - // forces - - ccelx = delx*ccel + fs1; - ccely = dely*ccel + fs2; - ccelz = delz*ccel + fs3; - - // stress contribution of atom pair to z-layers - // atom i always contributes - // atom j contributes if newton_pair is on or if owned by this proc - - m = static_cast ((x[i][dim]-boxlo) * stepinv); - if (m >= 0 && m < nlayers) { - sigxx[m] += delx*ccelx; - sigyy[m] += dely*ccely; - sigzz[m] += delz*ccelz; - sigxy[m] += delx*ccely; - sigxz[m] += delx*ccelz; - sigyz[m] += dely*ccelz; - } - - if (newton_pair || j < nlocal) { - m = static_cast ((x[j][dim]-boxlo) * stepinv); - if (m >= 0 && m < nlayers) { - sigxx[m] += delx*ccelx; - sigyy[m] += dely*ccely; - sigzz[m] += delz*ccelz; - sigxy[m] += delx*ccely; - sigxz[m] += delx*ccelz; - sigyz[m] += dely*ccelz; - } - } - } - } - } -} diff --git a/src/GRANULAR/fix_gran_diag.h b/src/GRANULAR/fix_gran_diag.h deleted file mode 100644 index 86fc8ec6b9..0000000000 --- a/src/GRANULAR/fix_gran_diag.h +++ /dev/null @@ -1,60 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef FIX_GRAN_DIAG_H -#define FIX_GRAN_DIAG_H - -#include "stdio.h" -#include "fix.h" - -namespace LAMMPS_NS { - -class FixGranDiag : public Fix { - public: - FixGranDiag(class LAMMPS *, int, char **); - ~FixGranDiag(); - int setmask(); - void init(); - void setup(); - void end_of_step(); - - private: - int me,first,pairstyle,nlayers,maxlayers,dim; - FILE *fpden,*fpvel,*fpstr; - double step,stepinv,boxlo,PI; - double dt,xkk,xkkt,xmu,gamman_dl,gammas_dl; - int freeze_group_bit; - - int *numdens; - double *dendens; - double *velx,*vely,*velz; - double *velxx,*velyy,*velzz,*velxy,*velxz,*velyz; - double *sigxx,*sigyy,*sigzz,*sigxy,*sigxz,*sigyz; - double *velx11,*vely11,*velz11; - double *velxx11,*velyy11,*velzz11,*velxy11,*velxz11,*velyz11; - double *velfxx,*velfyy,*velfzz,*velfxy,*velfxz,*velfyz; - double *sigx2,*sigy2,*sigz2,*sigxy2,*sigxz2,*sigyz2; - - class Pair *pair; - class NeighList *list,*listgran; - - void allocate(); - void deallocate(); - void stress_no_history(); - void stress_history(); - void stress_hertzian(); -}; - -} - -#endif diff --git a/src/GRANULAR/style_granular.h b/src/GRANULAR/style_granular.h index 7b0f4b6a71..5e89edb020 100644 --- a/src/GRANULAR/style_granular.h +++ b/src/GRANULAR/style_granular.h @@ -21,7 +21,6 @@ AtomStyle(granular,AtomVecGranular) #ifdef FixInclude #include "fix_freeze.h" -#include "fix_gran_diag.h" #include "fix_nve_gran.h" #include "fix_pour.h" #include "fix_shear_history.h" @@ -30,7 +29,6 @@ AtomStyle(granular,AtomVecGranular) #ifdef FixClass FixStyle(freeze,FixFreeze) -FixStyle(gran/diag,FixGranDiag) FixStyle(nve/gran,FixNVEGran) FixStyle(pour,FixPour) FixStyle(SHEAR_HISTORY,FixShearHistory) diff --git a/src/MEAM/pair_meam.cpp b/src/MEAM/pair_meam.cpp index d8031edf6e..9d52f2c988 100644 --- a/src/MEAM/pair_meam.cpp +++ b/src/MEAM/pair_meam.cpp @@ -122,6 +122,15 @@ void PairMEAM::compute(int eflag, int vflag) int *ilist_half,*jlist_half,*numneigh_half,**firstneigh_half; int *numneigh_full,**firstneigh_full; + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + // error check (for now) + + if (eflag_atom || vflag_atom) + error->all("Pair style meam does not yet support peratom energy/virial"); + // grow local arrays if necessary if (atom->nmax > nmax) { @@ -176,7 +185,7 @@ void PairMEAM::compute(int eflag, int vflag) // check size of scrfcn based on half neighbor list int nlocal = atom->nlocal; - int nall = atom->nlocal + atom->nghost; + int nall = nlocal + atom->nghost; n = 0; for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]]; @@ -194,9 +203,6 @@ void PairMEAM::compute(int eflag, int vflag) (double *) memory->smalloc(maxneigh*sizeof(double),"pair:fcpair"); } - eng_vdwl = 0.0; - if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; - // zero out local arrays for (i = 0; i < nall; i++) { @@ -229,7 +235,7 @@ void PairMEAM::compute(int eflag, int vflag) for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_dens_init_(&ifort,&nmax,&eflag,&eng_vdwl,&ntype,type,fmap,&x[0][0], + meam_dens_init_(&ifort,&nmax,&eflag,&evdwl,&ntype,type,fmap,&x[0][0], &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], @@ -246,7 +252,7 @@ void PairMEAM::compute(int eflag, int vflag) reverse_flag = 0; comm->reverse_comm_pair(this); - meam_dens_final_(&nlocal,&nmax,&eflag,&eng_vdwl,&ntype,type,fmap, + meam_dens_final_(&nlocal,&nmax,&eflag,&evdwl,&ntype,type,fmap, &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], &arho3b[0][0],&t_ave[0][0],gamma,dgamma1, dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); @@ -263,7 +269,7 @@ void PairMEAM::compute(int eflag, int vflag) for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_force_(&ifort,&nmax,&eflag,&eng_vdwl,&ntype,type,fmap,&x[0][0], + meam_force_(&ifort,&nmax,&eflag,&evdwl,&ntype,type,fmap,&x[0][0], &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], @@ -286,7 +292,10 @@ void PairMEAM::compute(int eflag, int vflag) neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half); neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full); - if (vflag == 2) virial_compute(); + // just sum global energy (for now) + + if (evflag) ev_tally(0,0,nlocal,newton_pair,evdwl,0.0,0.0,0.0,0.0); + if (vflag_fdotr) virial_compute(); } /* ---------------------------------------------------------------------- */ diff --git a/src/compute_ebond_atom.cpp b/src/compute_ebond_atom.cpp deleted file mode 100644 index 5cf7c845e1..0000000000 --- a/src/compute_ebond_atom.cpp +++ /dev/null @@ -1,166 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "string.h" -#include "compute_ebond_atom.h" -#include "atom.h" -#include "neighbor.h" -#include "comm.h" -#include "update.h" -#include "force.h" -#include "bond.h" -#include "modify.h" -#include "domain.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define MIN(a,b) ((a) < (b) ? (a) : (b)) -#define MAX(a,b) ((a) > (b) ? (a) : (b)) - -/* ---------------------------------------------------------------------- */ - -ComputeEbondAtom::ComputeEbondAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) -{ - if (narg != 3) error->all("Illegal compute ebond/atom command"); - - peratom_flag = 1; - size_peratom = 0; - comm_reverse = 1; - - nmax = 0; - energy = NULL; -} - -/* ---------------------------------------------------------------------- */ - -ComputeEbondAtom::~ComputeEbondAtom() -{ - memory->sfree(energy); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeEbondAtom::init() -{ - if (force->bond == NULL) - error->all("Bond style does not support computing per-atom bond energy"); - - int count = 0; - for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style,"ebond/atom") == 0) count++; - if (count > 1 && comm->me == 0) - error->warning("More than one compute ebond/atom"); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeEbondAtom::compute_peratom() -{ - int i,n,i1,i2,type,iflag,jflag; - double delx,dely,delz,rsq,fforce,eng; - - // grow energy array if necessary - - if (atom->nmax > nmax) { - memory->sfree(energy); - nmax = atom->nmax; - energy = (double *) - memory->smalloc(nmax*sizeof(double),"compute/epair/atom:energy"); - scalar_atom = energy; - } - - // clear energy array - // n includes ghosts only if newton_bond flag is set - - int nlocal = atom->nlocal; - int newton_bond = force->newton_bond; - - if (newton_bond) n = nlocal + atom->nghost; - else n = nlocal; - - for (i = 0; i < n; i++) energy[i] = 0.0; - - // compute bond energy for atoms via bond->single() - // if neither atom is in compute group, skip that bond - // only add energy to atoms in group - - double **x = atom->x; - int *mask = atom->mask; - int **bondlist = neighbor->bondlist; - int nbondlist = neighbor->nbondlist; - - for (n = 0; n < nbondlist; n++) { - i1 = bondlist[n][0]; - i2 = bondlist[n][1]; - type = bondlist[n][2]; - iflag = mask[i1] & groupbit; - jflag = mask[i2] & groupbit; - if (iflag == 0 && jflag == 0) continue; - - delx = x[i1][0] - x[i2][0]; - dely = x[i1][1] - x[i2][1]; - delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); - - rsq = delx*delx + dely*dely + delz*delz; - force->bond->single(type,rsq,i1,i2,1,fforce,eng); - if (iflag) energy[i1] += eng; - if (jflag && (newton_bond || i2 < nlocal)) energy[i2] += eng; - } - - // communicate energy between neigchbor procs - - if (newton_bond) comm->reverse_comm_compute(this); - - // remove double counting of per-atom energy - - for (i = 0; i < nlocal; i++) energy[i] *= 0.5; -} - -/* ---------------------------------------------------------------------- */ - -int ComputeEbondAtom::pack_reverse_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) buf[m++] = energy[i]; - return 1; -} - -/* ---------------------------------------------------------------------- */ - -void ComputeEbondAtom::unpack_reverse_comm(int n, int *list, double *buf) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - energy[j] += buf[m++]; - } -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based array -------------------------------------------------------------------------- */ - -double ComputeEbondAtom::memory_usage() -{ - double bytes = nmax * sizeof(double); - return bytes; -} diff --git a/src/compute_ebond_atom.h b/src/compute_ebond_atom.h deleted file mode 100644 index 6fb3388f71..0000000000 --- a/src/compute_ebond_atom.h +++ /dev/null @@ -1,38 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef COMPUTE_EBOND_ATOM_H -#define COMPUTE_EBOND_ATOM_H - -#include "compute.h" - -namespace LAMMPS_NS { - -class ComputeEbondAtom : public Compute { - public: - ComputeEbondAtom(class LAMMPS *, int, char **); - ~ComputeEbondAtom(); - void init(); - void compute_peratom(); - int pack_reverse_comm(int, int, double *); - void unpack_reverse_comm(int, int *, double *); - double memory_usage(); - - private: - int nmax; - double *energy; -}; - -} - -#endif diff --git a/src/compute_epair_atom.cpp b/src/compute_epair_atom.cpp deleted file mode 100644 index ae9b877684..0000000000 --- a/src/compute_epair_atom.cpp +++ /dev/null @@ -1,231 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "string.h" -#include "compute_epair_atom.h" -#include "atom.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "modify.h" -#include "comm.h" -#include "update.h" -#include "force.h" -#include "pair.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define MIN(a,b) ((a) < (b) ? (a) : (b)) -#define MAX(a,b) ((a) > (b) ? (a) : (b)) - -/* ---------------------------------------------------------------------- */ - -ComputeEpairAtom::ComputeEpairAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) -{ - if (narg != 3) error->all("Illegal compute epair/atom command"); - - peratom_flag = 1; - size_peratom = 0; - comm_reverse = 1; - - nmax = 0; - energy = NULL; -} - -/* ---------------------------------------------------------------------- */ - -ComputeEpairAtom::~ComputeEpairAtom() -{ - memory->sfree(energy); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeEpairAtom::init() -{ - if (force->pair == NULL || force->pair->single_enable == 0) - error->all("Pair style does not support computing per-atom energy"); - - // need an occasional half neighbor list - - int irequest = neighbor->request((void *) this); - neighbor->requests[irequest]->pair = 0; - neighbor->requests[irequest]->compute = 1; - neighbor->requests[irequest]->occasional = 1; - - if (force->pair_match("eam")) eamstyle = 1; - else eamstyle = 0; - - int count = 0; - for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style,"epair/atom") == 0) count++; - if (count > 1 && comm->me == 0) - error->warning("More than one compute epair/atom"); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeEpairAtom::init_list(int id, NeighList *ptr) -{ - list = ptr; -} - -/* ---------------------------------------------------------------------- */ - -void ComputeEpairAtom::compute_peratom() -{ - int i,j,ii,jj,n,inum,jnum,itype,jtype,iflag,jflag; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq; - double factor_coul,factor_lj,eng; - int *ilist,*jlist,*numneigh,**firstneigh; - Pair::One one; - - // grow energy array if necessary - - if (atom->nmax > nmax) { - memory->sfree(energy); - nmax = atom->nmax; - energy = (double *) - memory->smalloc(nmax*sizeof(double),"compute/epair/atom:energy"); - scalar_atom = energy; - } - - // clear energy array - // n includes ghosts only if newton_pair flag is set - - int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; - - if (newton_pair) n = nlocal + atom->nghost; - else n = nlocal; - - for (i = 0; i < n; i++) energy[i] = 0.0; - - // invoke half neighbor list (will copy or build if necessary) - - neighbor->build_one(list->index); - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // compute pairwise energy for atoms via pair->single() - // if neither atom is in compute group, skip that pair - // only add energy to atoms in group - - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - double **cutsq = force->pair->cutsq; - - double **x = atom->x; - int *mask = atom->mask; - int *type = atom->type; - int nall = nlocal + atom->nghost; - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - iflag = mask[i] & groupbit; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - jflag = mask[j] & groupbit; - if (iflag == 0 && jflag == 0) continue; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - force->pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,1,one); - eng = one.eng_coul + one.eng_vdwl; - if (iflag) energy[i] += eng; - if (jflag && (newton_pair || j < nlocal)) energy[j] += eng; - } - } - } - - // communicate energy between neigchbor procs - - if (newton_pair) comm->reverse_comm_compute(this); - - // remove double counting of per-atom energy - - for (i = 0; i < nlocal; i++) energy[i] *= 0.5; - - // for EAM, include embedding function contribution to energy - // only for atoms in compute group - - if (eamstyle) { - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - if (mask[i] & groupbit) { - force->pair->single_embed(i,type[i],eng); - energy[i] += eng; - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -int ComputeEpairAtom::pack_reverse_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) buf[m++] = energy[i]; - return 1; -} - -/* ---------------------------------------------------------------------- */ - -void ComputeEpairAtom::unpack_reverse_comm(int n, int *list, double *buf) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - energy[j] += buf[m++]; - } -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based array -------------------------------------------------------------------------- */ - -double ComputeEpairAtom::memory_usage() -{ - double bytes = nmax * sizeof(double); - return bytes; -} diff --git a/src/compute_epair_atom.h b/src/compute_epair_atom.h deleted file mode 100644 index 15c6667925..0000000000 --- a/src/compute_epair_atom.h +++ /dev/null @@ -1,40 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef COMPUTE_EPAIR_ATOM_H -#define COMPUTE_EPAIR_ATOM_H - -#include "compute.h" - -namespace LAMMPS_NS { - -class ComputeEpairAtom : public Compute { - public: - ComputeEpairAtom(class LAMMPS *, int, char **); - ~ComputeEpairAtom(); - void init(); - void init_list(int, class NeighList *); - void compute_peratom(); - int pack_reverse_comm(int, int, double *); - void unpack_reverse_comm(int, int *, double *); - double memory_usage(); - - private: - int nmax,eamstyle; - class NeighList *list; - double *energy; -}; - -} - -#endif