diff --git a/src/USER-MISC/compute_temp_rotate.cpp b/src/USER-MISC/compute_temp_rotate.cpp index 04e04a7520..02cc87c4b9 100644 --- a/src/USER-MISC/compute_temp_rotate.cpp +++ b/src/USER-MISC/compute_temp_rotate.cpp @@ -78,11 +78,9 @@ void ComputeTempRotate::setup() void ComputeTempRotate::dof_compute() { adjust_dof_fix(); - double natoms = group->count(igroup); - dof = domain->dimension * natoms; + natoms_temp = group->count(igroup); + dof = domain->dimension * natoms_temp; dof -= extra_dof + fix_dof; - if (dof < 0.0 && natoms > 0.0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -144,6 +142,8 @@ double ComputeTempRotate::compute_scalar() MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); + if (dof < 0.0 && natoms_temp > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } diff --git a/src/compute.h b/src/compute.h index ce0e7dc9e8..ba3835a34d 100644 --- a/src/compute.h +++ b/src/compute.h @@ -137,6 +137,7 @@ class Compute : protected Pointers { protected: int instance_me; // which Compute class instantiation I am + double natoms_temp; // # of atoms used for temperature calculation int extra_dof; // extra DOF for temperature computes int fix_dof; // DOF due to fixes int dynamic; // recount atoms for temperature computes diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp index 0e27aa4f5c..4a830d96af 100644 --- a/src/compute_temp.cpp +++ b/src/compute_temp.cpp @@ -61,11 +61,9 @@ void ComputeTemp::setup() void ComputeTemp::dof_compute() { adjust_dof_fix(); - double natoms = group->count(igroup); - dof = domain->dimension * natoms; + natoms_temp = group->count(igroup); + dof = domain->dimension * natoms_temp; dof -= extra_dof + fix_dof; - //if (dof < 0.0 && natoms > 0.0) - // error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -98,6 +96,8 @@ double ComputeTemp::compute_scalar() MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); + if (dof < 0.0 && natoms_temp > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } diff --git a/src/compute_temp_chunk.cpp b/src/compute_temp_chunk.cpp index afbe02129a..6713fcef85 100644 --- a/src/compute_temp_chunk.cpp +++ b/src/compute_temp_chunk.cpp @@ -288,10 +288,10 @@ double ComputeTempChunk::compute_scalar() MPI_Allreduce(&rcount,&allcount,1,MPI_DOUBLE,MPI_SUM,world); double dof = nchunk*cdof + adof*allcount; - if (dof < 0.0 && allcount > 0.0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); double tfactor = 0.0; if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); + if (dof < 0.0 && allcount > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp index e9f31dbfe1..9e734299af 100644 --- a/src/compute_temp_com.cpp +++ b/src/compute_temp_com.cpp @@ -70,11 +70,9 @@ void ComputeTempCOM::setup() void ComputeTempCOM::dof_compute() { adjust_dof_fix(); - double natoms = group->count(igroup); - dof = domain->dimension * natoms; + natoms_temp = group->count(igroup); + dof = domain->dimension * natoms_temp; dof -= extra_dof + fix_dof; - if (dof < 0.0 && natoms > 0.0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -113,6 +111,8 @@ double ComputeTempCOM::compute_scalar() MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); + if (dof < 0.0 && natoms_temp > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp index 7513f6fb1c..0035cbf4df 100644 --- a/src/compute_temp_deform.cpp +++ b/src/compute_temp_deform.cpp @@ -96,11 +96,9 @@ void ComputeTempDeform::setup() void ComputeTempDeform::dof_compute() { adjust_dof_fix(); - double natoms = group->count(igroup); - dof = domain->dimension * natoms; + natoms_temp = group->count(igroup); + dof = domain->dimension * natoms_temp; dof -= extra_dof + fix_dof; - if (dof < 0.0 && natoms > 0.0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -150,6 +148,8 @@ double ComputeTempDeform::compute_scalar() MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); + if (dof < 0.0 && natoms_temp > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp index 29525d3e80..f00a44f942 100644 --- a/src/compute_temp_partial.cpp +++ b/src/compute_temp_partial.cpp @@ -74,9 +74,9 @@ void ComputeTempPartial::setup() void ComputeTempPartial::dof_compute() { adjust_dof_fix(); - double natoms = group->count(igroup); + natoms_temp = group->count(igroup); int nper = xflag+yflag+zflag; - dof = nper * natoms; + dof = nper * natoms_temp; dof -= (1.0*nper/domain->dimension)*fix_dof + extra_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; @@ -119,7 +119,7 @@ double ComputeTempPartial::compute_scalar() MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); - if (tfactor == 0.0 && atom->natoms != 0) + if (dof < 0.0 && natoms_temp > 0.0) error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index 9f0ff3f31e..828c5c67a2 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -195,11 +195,9 @@ void ComputeTempProfile::setup() void ComputeTempProfile::dof_compute() { adjust_dof_fix(); - double natoms = group->count(igroup); - dof = domain->dimension * natoms; + natoms_temp = group->count(igroup); + dof = domain->dimension * natoms_temp; dof -= extra_dof + fix_dof; - if (dof < 0.0 && natoms > 0.0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -243,6 +241,8 @@ double ComputeTempProfile::compute_scalar() MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); + if (dof < 0.0 && natoms_temp > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp index 0c38ed811a..a922cded8c 100644 --- a/src/compute_temp_ramp.cpp +++ b/src/compute_temp_ramp.cpp @@ -126,11 +126,9 @@ void ComputeTempRamp::setup() void ComputeTempRamp::dof_compute() { adjust_dof_fix(); - double natoms = group->count(igroup); - dof = domain->dimension * natoms; + natoms_temp = group->count(igroup); + dof = domain->dimension * natoms_temp; dof -= extra_dof + fix_dof; - if (dof < 0.0 && natoms > 0.0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -172,6 +170,8 @@ double ComputeTempRamp::compute_scalar() MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); + if (dof < 0.0 && natoms_temp > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp index a90f7aefe2..e746ea1fc1 100644 --- a/src/compute_temp_region.cpp +++ b/src/compute_temp_region.cpp @@ -136,6 +136,8 @@ double ComputeTempRegion::compute_scalar() tarray[1] = t; MPI_Allreduce(tarray,tarray_all,2,MPI_DOUBLE,MPI_SUM,world); dof = domain->dimension * tarray_all[0] - extra_dof; + if (dof < 0.0 && tarray_all[0] > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) scalar = force->mvv2e * tarray_all[1] / (dof * force->boltz); else scalar = 0.0; return scalar; diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp index 050a987933..001a7c2de4 100644 --- a/src/compute_temp_sphere.cpp +++ b/src/compute_temp_sphere.cpp @@ -125,7 +125,7 @@ void ComputeTempSphere::dof_compute() int count,count_all; adjust_dof_fix(); - double natoms = group->count(igroup); + natoms_temp = group->count(igroup); // 6 or 3 dof for extended/point particles for 3d // 3 or 2 dof for extended/point particles for 2d @@ -166,7 +166,7 @@ void ComputeTempSphere::dof_compute() // additional adjustments to dof if (tempbias == 1) { - if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms; + if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms_temp; } else if (tempbias == 2) { int *mask = atom->mask; @@ -206,8 +206,6 @@ void ComputeTempSphere::dof_compute() } dof -= extra_dof + fix_dof; - if (dof < 0.0 && natoms > 0.0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -252,6 +250,8 @@ double ComputeTempSphere::compute_scalar() MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic || tempbias == 2) dof_compute(); + if (dof < 0.0 && natoms_temp > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; }