diff --git a/src/fix.h b/src/fix.h index 290ceebcda..af8e4dc554 100644 --- a/src/fix.h +++ b/src/fix.h @@ -47,7 +47,8 @@ class Fix : protected Pointers { int energy_atom_flag; // 1 if contributes to peratom eng int virial_global_flag; // 1 if contributes to global virial int virial_peratom_flag; // 1 if contributes to peratom virial - int ecouple_flag; // 1 if thermostat which accumulates eng to ecouple + int ecouple_flag; // 1 if thermostat fix outputs cumulative + // reservoir energy via compute_scalar() int time_integrate; // 1 if performs time integration, 0 if no int rigid_flag; // 1 if integrates rigid bodies, 0 if not int no_change_box; // 1 if cannot swap ortho <-> triclinic @@ -102,9 +103,8 @@ class Fix : protected Pointers { int comm_reverse; // size of reverse communication (0 if none) int comm_border; // size of border communication (0 if none) - double ecouple; // cumulative energy added to reservoir by thermostatting - double virial[6]; // virial for this timestep - double *eatom,**vatom; // per-atom energy/virial for this timestep + double virial[6]; // virial for this timestep + double *eatom,**vatom; // per-atom energy/virial for this timestep int centroidstressflag; // centroid stress compared to two-body stress // CENTROID_SAME = same as two-body stress diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp index f798252a50..cb7c87e6ff 100644 --- a/src/fix_langevin.cpp +++ b/src/fix_langevin.cpp @@ -53,7 +53,8 @@ enum{CONSTANT,EQUAL,ATOM}; FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), gjfflag(0), gfactor1(nullptr), gfactor2(nullptr), ratio(nullptr), tstr(nullptr), - flangevin(nullptr), tforce(nullptr), franprev(nullptr), lv(nullptr), id_temp(nullptr), random(nullptr) + flangevin(nullptr), tforce(nullptr), franprev(nullptr), + lv(nullptr), id_temp(nullptr), random(nullptr) { if (narg < 7) error->all(FLERR,"Illegal fix langevin command"); @@ -159,7 +160,7 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) : id_temp = nullptr; temperature = nullptr; - ecouple = 0.0; + energy = 0.0; // flangevin is unallocated until first call to setup() // compute_scalar checks for this and returns 0.0 @@ -990,7 +991,7 @@ void FixLangevin::end_of_step() } } - ecouple += energy_onestep*update->dt; + energy += energy_onestep*update->dt; } /* ---------------------------------------------------------------------- */ @@ -1068,7 +1069,7 @@ double FixLangevin::compute_scalar() if (mask[i] & groupbit) energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] + flangevin[i][2]*v[i][2]; - ecouple = 0.5*energy_onestep*update->dt; + energy = 0.5*energy_onestep*update->dt; } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -1079,13 +1080,13 @@ double FixLangevin::compute_scalar() if (tbiasflag) temperature->restore_bias(i, lv[i]); } - ecouple = -0.5*energy_onestep*update->dt; + energy = -0.5*energy_onestep*update->dt; } } // convert midstep energy back to previous fullstep energy - double energy_me = ecouple - 0.5*energy_onestep*update->dt; + double energy_me = energy - 0.5*energy_onestep*update->dt; double energy_all; MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index ee2b74ddb9..f4de07a5ae 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -51,7 +51,8 @@ enum{ISO,ANISO,TRICLINIC}; FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - rfix(nullptr), id_dilate(nullptr), irregular(nullptr), id_temp(nullptr), id_press(nullptr), + rfix(nullptr), id_dilate(nullptr), irregular(nullptr), + id_temp(nullptr), id_press(nullptr), eta(nullptr), eta_dot(nullptr), eta_dotdot(nullptr), eta_mass(nullptr), etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr), etap_mass(nullptr) @@ -784,7 +785,8 @@ void FixNH::setup(int /*vflag*/) } else { t0 = temperature->compute_scalar(); if (t0 < EPSILON) - error->all(FLERR, "Current temperature too close to zero, consider using ptemp setting"); + error->all(FLERR,"Current temperature too close to zero, " + "consider using ptemp setting"); } } t_target = t0; diff --git a/src/fix_temp_berendsen.cpp b/src/fix_temp_berendsen.cpp index 61cec81e8f..0ff5857fb5 100644 --- a/src/fix_temp_berendsen.cpp +++ b/src/fix_temp_berendsen.cpp @@ -81,7 +81,7 @@ FixTempBerendsen::FixTempBerendsen(LAMMPS *lmp, int narg, char **arg) : modify->add_compute(cmd); tflag = 1; - ecouple = 0; + energy = 0; } /* ---------------------------------------------------------------------- */ @@ -170,7 +170,7 @@ void FixTempBerendsen::end_of_step() double lamda = sqrt(1.0 + update->dt/t_period*(t_target/t_current - 1.0)); double efactor = 0.5 * force->boltz * tdof; - ecouple += t_current * (1.0-lamda*lamda) * efactor; + energy += t_current * (1.0-lamda*lamda) * efactor; double **v = atom->v; int *mask = atom->mask; @@ -238,7 +238,7 @@ void FixTempBerendsen::reset_target(double t_new) double FixTempBerendsen::compute_scalar() { - return ecouple; + return energy; } /* ---------------------------------------------------------------------- @@ -249,7 +249,7 @@ void FixTempBerendsen::write_restart(FILE *fp) { int n = 0; double list[1]; - list[n++] = ecouple; + list[n++] = energy; if (comm->me == 0) { int size = n * sizeof(double); @@ -266,7 +266,7 @@ void FixTempBerendsen::restart(char *buf) { double *list = (double *) buf; - ecouple = list[0]; + energy = list[0]; } /* ---------------------------------------------------------------------- diff --git a/src/fix_temp_csld.cpp b/src/fix_temp_csld.cpp index da6447b694..b4da730b97 100644 --- a/src/fix_temp_csld.cpp +++ b/src/fix_temp_csld.cpp @@ -92,7 +92,7 @@ FixTempCSLD::FixTempCSLD(LAMMPS *lmp, int narg, char **arg) : vhold = nullptr; nmax = -1; - ecouple = 0.0; + energy = 0.0; } /* ---------------------------------------------------------------------- */ @@ -250,7 +250,7 @@ void FixTempCSLD::end_of_step() // tally the kinetic energy transferred between heat bath and system t_current = temperature->compute_scalar(); - ecouple += ekin_old - t_current * 0.5 * temperature->dof * force->boltz; + energy += ekin_old - t_current * 0.5 * temperature->dof * force->boltz; } /* ---------------------------------------------------------------------- */ @@ -294,7 +294,7 @@ void FixTempCSLD::reset_target(double t_new) double FixTempCSLD::compute_scalar() { - return ecouple; + return energy; } /* ---------------------------------------------------------------------- @@ -304,11 +304,11 @@ double FixTempCSLD::compute_scalar() void FixTempCSLD::write_restart(FILE *fp) { const int PRNGSIZE = 98+2+3; - int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + ecouple + int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy double *list = nullptr; if (comm->me == 0) { list = new double[nsize]; - list[0] = ecouple; + list[0] = energy; list[1] = comm->nprocs; } double state[PRNGSIZE]; @@ -331,7 +331,7 @@ void FixTempCSLD::restart(char *buf) { double *list = (double *) buf; - ecouple = list[0]; + energy = list[0]; int nprocs = (int) list[1]; if (nprocs != comm->nprocs) { if (comm->me == 0) diff --git a/src/fix_temp_csvr.cpp b/src/fix_temp_csvr.cpp index 8779e7fb5e..b89cf9959a 100644 --- a/src/fix_temp_csvr.cpp +++ b/src/fix_temp_csvr.cpp @@ -92,7 +92,7 @@ FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) : tflag = 1; nmax = -1; - ecouple = 0.0; + energy = 0.0; } /* ---------------------------------------------------------------------- */ @@ -206,7 +206,7 @@ void FixTempCSVR::end_of_step() // tally the kinetic energy transferred between heat bath and system - ecouple += ekin_old * (1.0 - lamda*lamda); + energy += ekin_old * (1.0 - lamda*lamda); } /* ---------------------------------------------------------------------- */ @@ -329,7 +329,7 @@ void FixTempCSVR::reset_target(double t_new) double FixTempCSVR::compute_scalar() { - return ecouple; + return energy; } /* ---------------------------------------------------------------------- @@ -339,11 +339,11 @@ double FixTempCSVR::compute_scalar() void FixTempCSVR::write_restart(FILE *fp) { const int PRNGSIZE = 98+2+3; - int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + ecouple + int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy double *list = nullptr; if (comm->me == 0) { list = new double[nsize]; - list[0] = ecouple; + list[0] = energy; list[1] = comm->nprocs; } double state[PRNGSIZE]; @@ -366,7 +366,7 @@ void FixTempCSVR::restart(char *buf) { double *list = (double *) buf; - ecouple = list[0]; + energy = list[0]; int nprocs = (int) list[1]; if (nprocs != comm->nprocs) { if (comm->me == 0) diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp index 80470e3ccb..b393815faf 100644 --- a/src/fix_temp_rescale.cpp +++ b/src/fix_temp_rescale.cpp @@ -78,7 +78,7 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) : modify->add_compute(cmd); tflag = 1; - ecouple = 0.0; + energy = 0.0; } /* ---------------------------------------------------------------------- */ @@ -171,7 +171,7 @@ void FixTempRescale::end_of_step() int *mask = atom->mask; int nlocal = atom->nlocal; - ecouple += (t_current-t_target) * efactor; + energy += (t_current-t_target) * efactor; if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { @@ -236,7 +236,7 @@ void FixTempRescale::reset_target(double t_new) double FixTempRescale::compute_scalar() { - return ecouple; + return energy; } /* ---------------------------------------------------------------------- @@ -247,7 +247,7 @@ void FixTempRescale::write_restart(FILE *fp) { int n = 0; double list[1]; - list[n++] = ecouple; + list[n++] = energy; if (comm->me == 0) { int size = n * sizeof(double); @@ -265,7 +265,7 @@ void FixTempRescale::restart(char *buf) int n = 0; double *list = (double *) buf; - ecouple = list[n++]; + energy = list[n++]; } /* ---------------------------------------------------------------------- diff --git a/src/modify.cpp b/src/modify.cpp index 5812947dd4..f8718ca5df 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -489,7 +489,7 @@ void Modify::end_of_step() /* ---------------------------------------------------------------------- coupling energy call, only for relevant fixes - stored by each fix in ecouple variable + each thermostsat fix returns this via compute_scalar() ecouple = cumulative energy added to reservoir by thermostatting ------------------------------------------------------------------------- */ @@ -497,7 +497,7 @@ double Modify::energy_couple() { double energy = 0.0; for (int i = 0; i < n_energy_couple; i++) - energy += fix[list_energy_couple[i]]->ecouple; + energy += fix[list_energy_couple[i]]->compute_scalar(); return energy; }