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

This commit is contained in:
sjplimp 2013-01-09 18:14:19 +00:00
parent feefdd97ea
commit ec624ab446
2 changed files with 45 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

58
src/fix_heat.cpp
View File

@ -152,22 +152,22 @@ void FixHeat::end_of_step()
double *mass = atom->mass;
double *rmass = atom->rmass;
// reallocate vheat array if necessary
// reallocate per-atom arrays if necessary
if (hstyle == ATOM && atom->nlocal > maxatom) {
maxatom = atom->nmax;
memory->destroy(vheat);
memory->create(vheat,maxatom,"heat:vheat");
memory->destroy(vscale);
memory->create(vheat,maxatom,"heat:vheat");
memory->create(vscale,maxatom,"heat:vscale");
}
// evaluate variable
if (hstyle != CONSTANT) {
modify->clearstep_compute();
if (hstyle == EQUAL) heat_input = input->variable->compute_equal(hvar);
else input->variable->compute_atom(hvar,igroup,vheat,1,0);
modify->addstep_compute(update->ntimestep + nevery);
}
@ -184,19 +184,22 @@ void FixHeat::end_of_step()
}
double vcmsq = vcm[0]*vcm[0] + vcm[1]*vcm[1] + vcm[2]*vcm[2];
// scale = velocity scale factor to accomplish eflux change in energy
// vsub = velocity subtracted from each atom to preserve momentum
// add heat via scale factor on velocities for CONSTANT and EQUAL cases
// scale = velocity scale factor to accomplish eflux change in energy
// vsub = velocity subtracted from each atom to preserve momentum
// overall KE cannot go negative
if (hstyle != ATOM) {
heat = heat_input*nevery*update->dt*force->ftm2v;
double escale = (ke + heat - 0.5*vcmsq*masstotal)/(ke - 0.5*vcmsq*masstotal);
double escale =
(ke + heat - 0.5*vcmsq*masstotal)/(ke - 0.5*vcmsq*masstotal);
if (escale < 0.0) error->all(FLERR,"Fix heat kinetic energy went negative");
scale = sqrt(escale);
vsub[0] = (scale-1.0) * vcm[0];
vsub[1] = (scale-1.0) * vcm[1];
vsub[2] = (scale-1.0) * vcm[2];
if (iregion < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
@ -215,15 +218,21 @@ void FixHeat::end_of_step()
}
// add heat via per-atom scale factor on velocities for ATOM case
// vscale = velocity scale factor to accomplish eflux change in energy
// vsub = velocity subtracted from each atom to preserve momentum
// KE of an atom cannot go negative
} else {
vsub[0] = vsub[1] = vsub[2] = 0.0;
if (iregion < 0) {
for (i = 0; i < nlocal; i++)
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
heat = vheat[i]*nevery*update->dt*force->ftm2v;
vscale[i] = (ke + heat - 0.5*vcmsq*masstotal)/(ke - 0.5*vcmsq*masstotal);
if (vscale[i] < 0.0) error->all(FLERR,"Fix heat kinetic energy went negative");
vscale[i] =
(ke + heat - 0.5*vcmsq*masstotal)/(ke - 0.5*vcmsq*masstotal);
if (vscale[i] < 0.0)
error->all(FLERR,
"Fix heat kinetic energy of an atom went negative");
scale = sqrt(vscale[i]);
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
@ -231,9 +240,12 @@ void FixHeat::end_of_step()
vsub[1] += (scale-1.0) * v[i][1]*massone;
vsub[2] += (scale-1.0) * v[i][2]*massone;
}
}
vsub[0] /= masstotal;
vsub[1] /= masstotal;
vsub[2] /= masstotal;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
scale = sqrt(vscale[i]);
@ -241,13 +253,17 @@ void FixHeat::end_of_step()
v[i][1] = scale*v[i][1] - vsub[1];
v[i][2] = scale*v[i][2] - vsub[2];
}
} else {
region = domain->regions[iregion];
for (i = 0; i < nlocal; i++)
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
heat = vheat[i]*nevery*update->dt*force->ftm2v;
vscale[i] = (ke + heat - 0.5*vcmsq*masstotal)/(ke - 0.5*vcmsq*masstotal);
if (vscale[i] < 0.0) error->all(FLERR,"Fix heat kinetic energy went negative");
vscale[i] =
(ke + heat - 0.5*vcmsq*masstotal)/(ke - 0.5*vcmsq*masstotal);
if (vscale[i] < 0.0)
error->all(FLERR,
"Fix heat kinetic energy of an atom went negative");
scale = sqrt(vscale[i]);
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
@ -255,9 +271,12 @@ void FixHeat::end_of_step()
vsub[1] += (scale-1.0) * v[i][1]*massone;
vsub[2] += (scale-1.0) * v[i][2]*massone;
}
}
vsub[0] /= masstotal;
vsub[1] /= masstotal;
vsub[2] /= masstotal;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
scale = sqrt(vscale[i]);
@ -277,3 +296,14 @@ double FixHeat::compute_scalar()
return scale;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixHeat::memory_usage()
{
double bytes = 0.0;
if (hstyle == ATOM) bytes = atom->nmax*2 * sizeof(double);
return bytes;
}

1
src/fix_heat.h
View File

@ -32,6 +32,7 @@ class FixHeat : public Fix {
void init();
void end_of_step();
double compute_scalar();
double memory_usage();
private:
int iregion;