diff --git a/src/EXTRA-FIX/fix_ttm.cpp b/src/EXTRA-FIX/fix_ttm.cpp index 2af3685b4d..db7142e66b 100644 --- a/src/EXTRA-FIX/fix_ttm.cpp +++ b/src/EXTRA-FIX/fix_ttm.cpp @@ -287,16 +287,18 @@ void FixTTM::post_force(int /*vflag*/) int *mask = atom->mask; int nlocal = atom->nlocal; + double *boxlo = domain->boxlo; + double dxinv = nxgrid/domain->xprd; + double dyinv = nygrid/domain->yprd; + double dzinv = nzgrid/domain->zprd; + // apply damping and thermostat to all atoms in fix group for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - xscale = (x[i][0] - domain->boxlo[0])/domain->xprd; - yscale = (x[i][1] - domain->boxlo[1])/domain->yprd; - zscale = (x[i][2] - domain->boxlo[2])/domain->zprd; - ix = static_cast(xscale*nxgrid + shift) - OFFSET; - iy = static_cast(yscale*nygrid + shift) - OFFSET; - iz = static_cast(zscale*nzgrid + shift) - OFFSET; + ix = static_cast ((x[i][0]-boxlo[0])*dxinv + shift) - OFFSET; + iy = static_cast ((x[i][1]-boxlo[1])*dyinv + shift) - OFFSET; + iz = static_cast ((x[i][2]-boxlo[2])*dzinv + shift) - OFFSET; if (ix < 0) ix += nxgrid; if (iy < 0) iy += nygrid; if (iz < 0) iz += nzgrid; @@ -354,6 +356,11 @@ void FixTTM::end_of_step() int *mask = atom->mask; int nlocal = atom->nlocal; + double *boxlo = domain->boxlo; + double dxinv = nxgrid/domain->xprd; + double dyinv = nygrid/domain->yprd; + double dzinv = nzgrid/domain->zprd; + for (iz = 0; iz < nzgrid; iz++) for (iy = 0; iy < nygrid; iy++) for (ix = 0; ix < nxgrid; ix++) @@ -361,25 +368,19 @@ void FixTTM::end_of_step() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xscale = (x[i][0] - domain->boxlo[0])/domain->xprd; - yscale = (x[i][1] - domain->boxlo[1])/domain->yprd; - zscale = (x[i][2] - domain->boxlo[2])/domain->zprd; - ix = static_cast(xscale*nxgrid + shift) - OFFSET; - iy = static_cast(yscale*nygrid + shift) - OFFSET; - iz = static_cast(zscale*nzgrid + shift) - OFFSET; + ix = static_cast ((x[i][0]-boxlo[0])*dxinv + shift) - OFFSET; + iy = static_cast ((x[i][1]-boxlo[1])*dyinv + shift) - OFFSET; + iz = static_cast ((x[i][2]-boxlo[2])*dzinv + shift) - OFFSET; if (ix < 0) ix += nxgrid; if (iy < 0) iy += nygrid; if (iz < 0) iz += nzgrid; if (ix >= nxgrid) ix -= nxgrid; if (iy >= nygrid) iy -= nygrid; if (iz >= nzgrid) iz -= nzgrid; + net_energy_transfer[iz][iy][ix] += (flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] + flangevin[i][2]*v[i][2]); - //printf("FLANG i %d xyz %g %g %g\n",i, - // flangevin[i][0], - // flangevin[i][1], - // flangevin[i][2]); } outflag = 0; @@ -443,15 +444,12 @@ void FixTTM::end_of_step() inner_dt/(electronic_specific_heat*electronic_density) * (electronic_thermal_conductivity * - ((T_electron_old[iz][iy][xright] + - T_electron_old[iz][iy][xleft] - - 2*T_electron_old[iz][iy][ix])/dx/dx + - (T_electron_old[iz][yright][ix] + - T_electron_old[iz][yleft][ix] - - 2*T_electron_old[iz][iy][ix])/dy/dy + - (T_electron_old[zright][iy][ix] + - T_electron_old[zleft][iy][ix] - - 2*T_electron_old[iz][iy][ix])/dz/dz) - + ((T_electron_old[iz][iy][xright] + T_electron_old[iz][iy][xleft] - + 2.0*T_electron_old[iz][iy][ix])/dx/dx + + (T_electron_old[iz][yright][ix] + T_electron_old[iz][yleft][ix] - + 2.0*T_electron_old[iz][iy][ix])/dy/dy + + (T_electron_old[zright][iy][ix] + T_electron_old[zleft][iy][ix] - + 2.0*T_electron_old[iz][iy][ix])/dz/dz) - (net_energy_transfer_all[iz][iy][ix])/del_vol); } diff --git a/src/EXTRA-FIX/fix_ttm_grid.cpp b/src/EXTRA-FIX/fix_ttm_grid.cpp index 4ded3448d6..4e7b06c668 100644 --- a/src/EXTRA-FIX/fix_ttm_grid.cpp +++ b/src/EXTRA-FIX/fix_ttm_grid.cpp @@ -79,8 +79,13 @@ void FixTTMGrid::post_constructor() ngridout*sizeof(double)); // set initial electron temperatures from user input file + // communicate new T_electron values to ghost grid points - if (infile) read_electron_temperatures(infile); + if (infile) { + read_electron_temperatures(infile); + gc->forward_comm(GridComm::FIX,this,1,sizeof(double),0, + gc_buf1,gc_buf2,MPI_DOUBLE); + } } /* ---------------------------------------------------------------------- */ @@ -118,7 +123,6 @@ void FixTTMGrid::post_force(int /*vflag*/) for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - ix = static_cast ((x[i][0]-boxlo[0])*dxinv + shift) - OFFSET; iy = static_cast ((x[i][1]-boxlo[1])*dyinv + shift) - OFFSET; iz = static_cast ((x[i][2]-boxlo[2])*dzinv + shift) - OFFSET; @@ -183,6 +187,7 @@ void FixTTMGrid::end_of_step() ix = static_cast ((x[i][0]-boxlo[0])*dxinv + shift) - OFFSET; iy = static_cast ((x[i][1]-boxlo[1])*dyinv + shift) - OFFSET; iz = static_cast ((x[i][2]-boxlo[2])*dzinv + shift) - OFFSET; + net_energy_transfer[iz][iy][ix] += (flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] + flangevin[i][2]*v[i][2]); @@ -228,21 +233,21 @@ void FixTTMGrid::end_of_step() (electronic_thermal_conductivity * ((T_electron_old[iz][iy][ix-1] + T_electron_old[iz][iy][ix+1] - - 2*T_electron_old[iz][iy][ix])*dxinv*dxinv + + 2.0*T_electron_old[iz][iy][ix])*dxinv*dxinv + (T_electron_old[iz][iy-1][ix] + T_electron_old[iz][iy+1][ix] - - 2*T_electron_old[iz][iy][ix])*dyinv*dyinv + + 2.0*T_electron_old[iz][iy][ix])*dyinv*dyinv + (T_electron_old[iz-1][iy][ix] + T_electron_old[iz+1][iy][ix] - - 2*T_electron_old[iz][iy][ix])*dzinv*dzinv) - + 2.0*T_electron_old[iz][iy][ix])*dzinv*dzinv) - net_energy_transfer[iz][iy][ix]/volgrid); } + + // communicate new T_electron values to ghost grid points + + gc->forward_comm(GridComm::FIX,this,1,sizeof(double),0, + gc_buf1,gc_buf2,MPI_DOUBLE); } - // communicate new T_electron values to ghost grid points - - gc->forward_comm(GridComm::FIX,this,1,sizeof(double),0, - gc_buf1,gc_buf2,MPI_DOUBLE); - // output of grid temperatures to file if (outfile && (update->ntimestep % outevery == 0)) { @@ -364,11 +369,6 @@ void FixTTMGrid::read_electron_temperatures(const char *filename) error->all(FLERR,"Fix ttm/grid infile did not set all temperatures"); memory->destroy3d_offset(T_initial_set,nzlo_in,nylo_in,nxlo_in); - - // communicate new T_electron values to ghost grid points - - gc->forward_comm(GridComm::FIX,this,1,sizeof(double),0, - gc_buf1,gc_buf2,MPI_DOUBLE); } /* ---------------------------------------------------------------------- @@ -593,6 +593,9 @@ void FixTTMGrid::restart(char *buf) for (ix = nxlo_in; ix <= nxhi_in; ix++) { iglobal = nygrid*nxgrid*iz + nxgrid*iy + ix; T_electron[iz][iy][ix] = rlist[n+iglobal]; + if (ix == 0 && iy == 0 & iz == 0) + printf("READRESTART 000 n+iglobal %d %d value %g\n",n,iglobal, + T_electron[ix][iy][iz]); } }