diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index 61ba36ea2b..cb2daaac2a 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -33,8 +33,6 @@ #include "memory.h" #include "error.h" - - using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; diff --git a/src/compute_pe.cpp b/src/compute_pe.cpp index f10f1eb344..b7078675e5 100644 --- a/src/compute_pe.cpp +++ b/src/compute_pe.cpp @@ -99,7 +99,8 @@ double ComputePE::compute_scalar() scalar += force->pair->etail / volume; } - if (fixflag && modify->n_thermo_energy) scalar += modify->thermo_energy(); + if (fixflag && modify->n_energy_global) + scalar += modify->energy_global(); return scalar; } diff --git a/src/compute_pe_atom.cpp b/src/compute_pe_atom.cpp index 52e132a626..ad490c66a6 100644 --- a/src/compute_pe_atom.cpp +++ b/src/compute_pe_atom.cpp @@ -152,8 +152,8 @@ void ComputePEAtom::compute_peratom() // add in per-atom contributions from relevant fixes // always only for owned atoms, not ghost - if (fixflag && modify->n_thermo_energy_atom) - modify->thermo_energy_atom(nlocal,energy); + if (fixflag && modify->n_energy_atom) + modify->energy_atom(nlocal,energy); // communicate ghost energy between neighbor procs diff --git a/src/fix.cpp b/src/fix.cpp index 573d8ae1a1..bd3d91767b 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -63,9 +63,11 @@ Fix::Fix(LAMMPS *lmp, int /*narg*/, char **arg) : box_change = NO_BOX_CHANGE; thermo_energy = 0; thermo_virial = 0; + energy_global_flag = energy_atom_flag = 0; + virial_global_flag = virial_atom_flag = 0; + ecouple_flag = 0; rigid_flag = 0; peatom_flag = 0; - virial_flag = 0; no_change_box = 0; time_integrate = 0; time_depend = 0; @@ -150,8 +152,7 @@ void Fix::modify_params(int narg, char **arg) if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command"); if (strcmp(arg[iarg+1],"no") == 0) thermo_energy = 0; else if (strcmp(arg[iarg+1],"yes") == 0) { - if (!(THERMO_ENERGY & setmask())) - error->all(FLERR,"Illegal fix_modify command"); + if (energy_flag == 0) error->all(FLERR,"Illegal fix_modify command"); thermo_energy = 1; } else error->all(FLERR,"Illegal fix_modify command"); iarg += 2; @@ -159,8 +160,7 @@ void Fix::modify_params(int narg, char **arg) if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command"); if (strcmp(arg[iarg+1],"no") == 0) thermo_virial = 0; else if (strcmp(arg[iarg+1],"yes") == 0) { - if (virial_flag == 0) - error->all(FLERR,"Illegal fix_modify command"); + if (virial_flag == 0) error->all(FLERR,"Illegal fix_modify command"); thermo_virial = 1; } else error->all(FLERR,"Illegal fix_modify command"); iarg += 2; @@ -180,9 +180,12 @@ void Fix::modify_params(int narg, char **arg) } /* ---------------------------------------------------------------------- - setup for energy, virial computation + setup for peratom energy and global/peratom virial computation see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) - fixes call this if they use ev_tally() + fixes call Fix::ev_init() if tally energy and virial values + if thermo_energy is not set, energy tallying is disabled + if thermo_virial is not set, virial tallying is disabled + global energy is tallied separately, output by compute_scalar() method ------------------------------------------------------------------------- */ void Fix::ev_setup(int eflag, int vflag) @@ -191,13 +194,19 @@ void Fix::ev_setup(int eflag, int vflag) evflag = 1; - eflag_either = eflag; - eflag_global = eflag & ENERGY_GLOBAL; - eflag_atom = eflag & ENERGY_ATOM; + if (!thermo_energy) eflag_either = eflag_global = eflag_atom = 0; + else { + eflag_either = eflag; + eflag_global = eflag & ENERGY_GLOBAL; + eflag_atom = eflag & ENERGY_ATOM; + } - vflag_either = vflag; - vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); - vflag_atom = vflag & (VIRIAL_ATOM | VIRIAL_CENTROID); + if (!thermo_virial) vflag_either = vflag_global = vflag_atom = 0; + else { + vflag_either = vflag; + vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); + vflag_atom = vflag & (VIRIAL_ATOM | VIRIAL_CENTROID); + } // reallocate per-atom arrays if necessary @@ -235,24 +244,17 @@ void Fix::ev_setup(int eflag, int vflag) } /* ---------------------------------------------------------------------- - if thermo_virial is on: - setup for virial computation - see integrate::ev_set() for values of vflag - fixes call this if use v_tally() - else: set evflag=0 + setup for global/peratom virial computation + see integrate::ev_set() for values of vflag (0-6) + fixes call Fix::v_init() if tally virial values but not energy + if thermo_virial is not set, virial tallying is disabled ------------------------------------------------------------------------- */ void Fix::v_setup(int vflag) { int i,n; - if (!thermo_virial) { - evflag = 0; - return; - } - evflag = 1; - vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); vflag_atom = vflag & (VIRIAL_ATOM | VIRIAL_CENTROID); @@ -304,7 +306,6 @@ void Fix::ev_tally(int n, int *list, double total, double eng, double *v) v_tally(n,list,total,v); } - /* ---------------------------------------------------------------------- tally virial into global and per-atom accumulators n = # of local owned atoms involved, with local indices in list diff --git a/src/fix.h b/src/fix.h index ec951e3bb2..290ceebcda 100644 --- a/src/fix.h +++ b/src/fix.h @@ -40,14 +40,17 @@ class Fix : protected Pointers { }; bigint next_reneighbor; // next timestep to force a reneighboring - int thermo_energy; // 1 if fix_modify enabled ThEng, 0 if not - int thermo_virial; // 1 if fix_modify enabled ThVir, 0 if not int nevery; // how often to call an end_of_step fix - int rigid_flag; // 1 if Fix integrates rigid bodies, 0 if not - int peatom_flag; // 1 if Fix contributes per-atom eng, 0 if not - int virial_flag; // 1 if Fix contributes to virial, 0 if not + int thermo_energy; // 1 if fix_modify energy enabled, 0 if not + int thermo_virial; // 1 if fix_modify virial enabled, 0 if not + int energy_global_flag; // 1 if contributes to global eng + 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 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 - int time_integrate; // 1 if fix performs time integration, 0 if no int time_depend; // 1 if requires continuous timestepping int create_attribute; // 1 if fix stores attributes that need // setting when a new atom is created @@ -99,8 +102,9 @@ 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 virial[6]; // accumulated virial - double *eatom,**vatom; // accumulated per-atom energy/virial + 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 int centroidstressflag; // centroid stress compared to two-body stress // CENTROID_SAME = same as two-body stress @@ -239,11 +243,17 @@ class Fix : protected Pointers { int dynamic; // recount atoms for temperature computes void ev_init(int eflag, int vflag) { - if (eflag||vflag) ev_setup(eflag, vflag); - else evflag = eflag_either = eflag_global = eflag_atom = vflag_either = vflag_global = vflag_atom = 0; + if ((eflag && thermo_energy) || (vflag && thermo_virial)) ev_setup(eflag, vflag); + else evflag = eflag_either = eflag_global = eflag_atom = + vflag_either = vflag_global = vflag_atom = 0; } void ev_setup(int, int); void ev_tally(int, int *, double, double, double *); + + void v_init(int vflag) { + if (vflag && thermo_virial) v_setup(vflag); + else evflag = vflag_either = vflag_global = vflag_atom = 0; + } void v_setup(int); void v_tally(int, int *, double, double *); void v_tally(int, double *); @@ -263,7 +273,7 @@ namespace FixConst { FINAL_INTEGRATE = 1<<8, END_OF_STEP = 1<<9, POST_RUN = 1<<10, - THERMO_ENERGY = 1<<11, + //THERMO_ENERGY = 1<<11, // remove when done with refactoring INITIAL_INTEGRATE_RESPA = 1<<12, POST_INTEGRATE_RESPA = 1<<13, PRE_FORCE_RESPA = 1<<14, diff --git a/src/fix_addforce.cpp b/src/fix_addforce.cpp index 881f77c05c..8bd5bc69a3 100644 --- a/src/fix_addforce.cpp +++ b/src/fix_addforce.cpp @@ -352,7 +352,7 @@ void FixAddForce::post_force(int vflag) v[3] = xstyle ? xvalue*unwrap[1] : 0.0; v[4] = xstyle ? xvalue*unwrap[2] : 0.0; v[5] = ystyle ? yvalue*unwrap[2] : 0.0; - v_tally(i, v); + v_tally(i,v); } } } diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp index 46b00f9e83..5e6638a56c 100644 --- a/src/fix_temp_rescale.cpp +++ b/src/fix_temp_rescale.cpp @@ -12,9 +12,10 @@ ------------------------------------------------------------------------- */ #include "fix_temp_rescale.h" -#include +#include #include + #include "atom.h" #include "force.h" #include "group.h" @@ -26,7 +27,6 @@ #include "compute.h" #include "error.h" - using namespace LAMMPS_NS; using namespace FixConst; diff --git a/src/modify.cpp b/src/modify.cpp index f3ebb03c38..5812947dd4 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -46,8 +46,8 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) n_initial_integrate = n_post_integrate = 0; n_pre_exchange = n_pre_neighbor = n_post_neighbor = 0; n_pre_force = n_pre_reverse = n_post_force = 0; - n_final_integrate = n_end_of_step = n_thermo_energy = 0; - n_thermo_energy_atom = 0; + n_final_integrate = n_end_of_step = 0; + n_energy_couple = n_energy_global = n_energy_atom = 0; n_initial_integrate_respa = n_post_integrate_respa = 0; n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0; n_min_pre_exchange = n_min_pre_force = n_min_pre_reverse = 0; @@ -60,7 +60,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) list_pre_exchange = list_pre_neighbor = list_post_neighbor = nullptr; list_pre_force = list_pre_reverse = list_post_force = nullptr; list_final_integrate = list_end_of_step = nullptr; - list_thermo_energy = list_thermo_energy_atom = nullptr; + list_energy_couple = list_energy_global = list_energy_atom = nullptr; list_initial_integrate_respa = list_post_integrate_respa = nullptr; list_pre_force_respa = list_post_force_respa = nullptr; list_final_integrate_respa = nullptr; @@ -137,8 +137,9 @@ Modify::~Modify() delete [] list_post_force; delete [] list_final_integrate; delete [] list_end_of_step; - delete [] list_thermo_energy; - delete [] list_thermo_energy_atom; + delete [] list_energy_couple; + delete [] list_energy_global; + delete [] list_energy_atom; delete [] list_initial_integrate_respa; delete [] list_post_integrate_respa; delete [] list_pre_force_respa; @@ -218,8 +219,9 @@ void Modify::init() list_init(POST_FORCE,n_post_force,list_post_force); list_init(FINAL_INTEGRATE,n_final_integrate,list_final_integrate); list_init_end_of_step(END_OF_STEP,n_end_of_step,list_end_of_step); - list_init_thermo_energy(THERMO_ENERGY,n_thermo_energy,list_thermo_energy); - list_init_thermo_energy_atom(n_thermo_energy_atom,list_thermo_energy_atom); + list_init_energy_couple(n_energy_couple,list_energy_couple); + list_init_energy_global(n_energy_global,list_energy_global); + list_init_energy_atom(n_energy_atom,list_energy_atom); list_init(INITIAL_INTEGRATE_RESPA, n_initial_integrate_respa,list_initial_integrate_respa); @@ -486,31 +488,45 @@ void Modify::end_of_step() } /* ---------------------------------------------------------------------- - thermo energy call, only for relevant fixes - called by Thermo class - compute_scalar() is fix call to return energy + coupling energy call, only for relevant fixes + stored by each fix in ecouple variable + ecouple = cumulative energy added to reservoir by thermostatting ------------------------------------------------------------------------- */ -double Modify::thermo_energy() +double Modify::energy_couple() { double energy = 0.0; - for (int i = 0; i < n_thermo_energy; i++) - energy += fix[list_thermo_energy[i]]->compute_scalar(); + for (int i = 0; i < n_energy_couple; i++) + energy += fix[list_energy_couple[i]]->ecouple; return energy; } /* ---------------------------------------------------------------------- - per-atom thermo energy call, only for relevant fixes + global energy call, only for relevant fixes + they return energy via compute_scalar() + called by compute pe +------------------------------------------------------------------------- */ + +double Modify::energy_global() +{ + double energy = 0.0; + for (i = 0; i < n_energy_global; i++) + energy += fix[list_energy_global[i]]->compute_scalar(); + return energy; +} + +/* ---------------------------------------------------------------------- + peratom energy call, only for relevant fixes called by compute pe/atom ------------------------------------------------------------------------- */ -void Modify::thermo_energy_atom(int nlocal, double *energy) +void Modify::energy_atom(int nlocal, double *energy) { int i,j; double *eatom; - for (i = 0; i < n_thermo_energy_atom; i++) { - eatom = fix[list_thermo_energy_atom[i]]->eatom; + for (i = 0; i < n_energy_atom; i++) { + eatom = fix[list_energy_atom[i]]->eatom; if (!eatom) continue; for (j = 0; j < nlocal; j++) energy[j] += eatom[j]; } @@ -1632,43 +1648,79 @@ void Modify::list_init_end_of_step(int mask, int &n, int *&list) } /* ---------------------------------------------------------------------- - create list of fix indices for thermo energy fixes - only added to list if fix has THERMO_ENERGY mask set, - and its thermo_energy flag was set via fix_modify + create list of fix indices for fixes that compute reservoir coupling energy + only added to list if fix has ecouple_flag set ------------------------------------------------------------------------- */ -void Modify::list_init_thermo_energy(int mask, int &n, int *&list) +void Modify::list_init_energy_couple(int &n, int *&list) { delete [] list; n = 0; for (int i = 0; i < nfix; i++) - if (fmask[i] & mask && fix[i]->thermo_energy) n++; + if (fix[i]->ecouple_flag) n++; list = new int[n]; n = 0; for (int i = 0; i < nfix; i++) - if (fmask[i] & mask && fix[i]->thermo_energy) list[n++] = i; + if (fix[i]->ecouple_flag) list[n++] = i; } /* ---------------------------------------------------------------------- - create list of fix indices for peratom thermo energy fixes - only added to list if fix has its peatom_flag set, - and its thermo_energy flag was set via fix_modify + create list of fix indices for fixes that compute peratom energy + only added to list if fix has energy_atom_flag and thermo_energy set ------------------------------------------------------------------------- */ -void Modify::list_init_thermo_energy_atom(int &n, int *&list) +void Modify::list_init_energy_atom(int &n, int *&list) { delete [] list; n = 0; for (int i = 0; i < nfix; i++) - if (fix[i]->peatom_flag && fix[i]->thermo_energy) n++; + if (fix[i]->energy_atom_flag && fix[i]->thermo_energy) n++; list = new int[n]; n = 0; for (int i = 0; i < nfix; i++) - if (fix[i]->peatom_flag && fix[i]->thermo_energy) list[n++] = i; + if (fix[i]->energy_atom_flag && fix[i]->thermo_energy) list[n++] = i; +} + +/* ---------------------------------------------------------------------- + create list of fix indices for fixes that compute global energy + only added to list if fix has energy_global_flag and thermo_energy set +------------------------------------------------------------------------- */ + +void Modify::list_init_energy_global(int &n, int *&list) +{ + delete [] list; + + n = 0; + for (int i = 0; i < nfix; i++) + if (fix[i]->energy_global_flag && fix[i]->thermo_energy) n++; + list = new int[n]; + + n = 0; + for (int i = 0; i < nfix; i++) + if (fix[i]->energy_global_flag && fix[i]->thermo_energy) list[n++] = i; +} + +/* ---------------------------------------------------------------------- + create list of fix indices for fixes that compute peratom energy + only added to list if fix has energy_atom_flag and thermo_energy set +------------------------------------------------------------------------- */ + +void Modify::list_init_energy_atom(int &n, int *&list) +{ + delete [] list; + + n = 0; + for (int i = 0; i < nfix; i++) + if (fix[i]->energy_atom_flag && fix[i]->thermo_energy) n++; + list = new int[n]; + + n = 0; + for (int i = 0; i < nfix; i++) + if (fix[i]->energy_atom_flag && fix[i]->thermo_energy) list[n++] = i; } /* ---------------------------------------------------------------------- diff --git a/src/modify.h b/src/modify.h index a347e8486d..ba8efd6525 100644 --- a/src/modify.h +++ b/src/modify.h @@ -33,7 +33,8 @@ class Modify : protected Pointers { int n_initial_integrate,n_post_integrate,n_pre_exchange; int n_pre_neighbor,n_post_neighbor; int n_pre_force,n_pre_reverse,n_post_force; - int n_final_integrate,n_end_of_step,n_thermo_energy,n_thermo_energy_atom; + int n_final_integrate,n_end_of_step; + int n_energy_couple,n_energy_global,n_energy_atom; int n_initial_integrate_respa,n_post_integrate_respa; int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa; int n_min_pre_exchange,n_min_pre_neighbor,n_min_post_neighbor; @@ -69,8 +70,9 @@ class Modify : protected Pointers { virtual void post_force(int); virtual void final_integrate(); virtual void end_of_step(); - virtual double thermo_energy(); - virtual void thermo_energy_atom(int, double *); + virtual double energy_couple(); + virtual double energy_global(); + virtual void energy_atom(int, double *); virtual void post_run(); virtual void create_attribute(int); @@ -135,8 +137,8 @@ class Modify : protected Pointers { int *list_initial_integrate,*list_post_integrate; int *list_pre_exchange,*list_pre_neighbor,*list_post_neighbor; int *list_pre_force,*list_pre_reverse,*list_post_force; - int *list_final_integrate,*list_end_of_step,*list_thermo_energy; - int *list_thermo_energy_atom; + int *list_final_integrate,*list_end_of_step; + int *list_energy_couple,*list_energy_global,*list_energy_atom; int *list_initial_integrate_respa,*list_post_integrate_respa; int *list_pre_force_respa,*list_post_force_respa; int *list_final_integrate_respa; @@ -163,8 +165,9 @@ class Modify : protected Pointers { void list_init(int, int &, int *&); void list_init_end_of_step(int, int &, int *&); - void list_init_thermo_energy(int, int &, int *&); - void list_init_thermo_energy_atom(int &, int *&); + void list_init_energy_couple(int &, int *&); + void list_init_energy_global(int &, int *&); + void list_init_energy_atom(int &, int *&); void list_init_dofflag(int &, int *&); void list_init_compute(); diff --git a/src/thermo.cpp b/src/thermo.cpp index 9197f88084..626dc09c3f 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -54,13 +54,14 @@ using namespace MathConst; // step, elapsed, elaplong, dt, time, cpu, tpcpu, spcpu, cpuremain, // part, timeremain -// atoms, temp, press, pe, ke, etotal, enthalpy +// atoms, temp, press, pe, ke, etotal // evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail -// vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz, +// enthalpy, ecouple, econserve +// vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz // xlat, ylat, zlat -// bonds, angles, dihedrals, impropers, +// bonds, angles, dihedrals, impropers // pxx, pyy, pzz, pxy, pxz, pyz -// fmax, fnorm, nbuild, ndanger, +// fmax, fnorm, nbuild, ndanger // cella, cellb, cellc, cellalpha, cellbeta, cellgamma // customize a new thermo style by adding a DEFINE to this list @@ -753,11 +754,6 @@ void Thermo::parse_fields(char *str) addfield("TotEng",&Thermo::compute_etotal,FLOAT); index_temp = add_compute(id_temp,SCALAR); index_pe = add_compute(id_pe,SCALAR); - } else if (word == "enthalpy") { - addfield("Enthalpy",&Thermo::compute_enthalpy,FLOAT); - index_temp = add_compute(id_temp,SCALAR); - index_press_scalar = add_compute(id_press,SCALAR); - index_pe = add_compute(id_pe,SCALAR); } else if (word == "evdwl") { addfield("E_vdwl",&Thermo::compute_evdwl,FLOAT); @@ -790,6 +786,19 @@ void Thermo::parse_fields(char *str) addfield("E_tail",&Thermo::compute_etail,FLOAT); index_pe = add_compute(id_pe,SCALAR); + } else if (word == "enthalpy") { + addfield("Enthalpy",&Thermo::compute_enthalpy,FLOAT); + index_temp = add_compute(id_temp,SCALAR); + index_press_scalar = add_compute(id_press,SCALAR); + index_pe = add_compute(id_pe,SCALAR); + } else if (word == "ecouple") { + addfield("Ecouple",&Thermo::compute_ecouple,FLOAT); + index_pe = add_compute(id_pe,SCALAR); + } else if (word == "econserve") { + addfield("Econserve",&Thermo::compute_econserve,FLOAT); + index_temp = add_compute(id_temp,SCALAR); + index_pe = add_compute(id_pe,SCALAR); + } else if (word == "vol") { addfield("Volume",&Thermo::compute_vol,FLOAT); } else if (word == "density") { @@ -1234,42 +1243,6 @@ int Thermo::evaluate_keyword(const char *word, double *answer) } compute_etotal(); - } else if (strcmp(word,"enthalpy") == 0) { - if (!pe) - error->all(FLERR, - "Thermo keyword in variable requires thermo to use/init pe"); - if (update->whichflag == 0) { - if (pe->invoked_scalar != update->ntimestep) - error->all(FLERR,"Compute used in variable thermo keyword between runs " - "is not current"); - } else { - pe->compute_scalar(); - pe->invoked_flag |= INVOKED_SCALAR; - } - if (!temperature) - error->all(FLERR,"Thermo keyword in variable requires " - "thermo to use/init temp"); - if (update->whichflag == 0) { - if (temperature->invoked_scalar != update->ntimestep) - error->all(FLERR,"Compute used in variable thermo keyword between runs " - "is not current"); - } else if (!(temperature->invoked_flag & INVOKED_SCALAR)) { - temperature->compute_scalar(); - temperature->invoked_flag |= INVOKED_SCALAR; - } - if (!pressure) - error->all(FLERR,"Thermo keyword in variable requires " - "thermo to use/init press"); - if (update->whichflag == 0) { - if (pressure->invoked_scalar != update->ntimestep) - error->all(FLERR,"Compute used in variable thermo keyword between runs " - "is not current"); - } else if (!(pressure->invoked_flag & INVOKED_SCALAR)) { - pressure->compute_scalar(); - pressure->invoked_flag |= INVOKED_SCALAR; - } - compute_enthalpy(); - } else if (strcmp(word,"evdwl") == 0) { if (update->eflag_global != update->ntimestep) error->all(FLERR,"Energy was not tallied on needed timestep"); @@ -1356,6 +1329,68 @@ int Thermo::evaluate_keyword(const char *word, double *answer) error->all(FLERR,"Energy was not tallied on needed timestep"); compute_etail(); + } else if (strcmp(word,"enthalpy") == 0) { + if (!pe) + error->all(FLERR, + "Thermo keyword in variable requires thermo to use/init pe"); + if (update->whichflag == 0) { + if (pe->invoked_scalar != update->ntimestep) + error->all(FLERR,"Compute used in variable thermo keyword between runs " + "is not current"); + } else { + pe->compute_scalar(); + pe->invoked_flag |= INVOKED_SCALAR; + } + if (!temperature) + error->all(FLERR,"Thermo keyword in variable requires " + "thermo to use/init temp"); + if (update->whichflag == 0) { + if (temperature->invoked_scalar != update->ntimestep) + error->all(FLERR,"Compute used in variable thermo keyword between runs " + "is not current"); + } else if (!(temperature->invoked_flag & INVOKED_SCALAR)) { + temperature->compute_scalar(); + temperature->invoked_flag |= INVOKED_SCALAR; + } + if (!pressure) + error->all(FLERR,"Thermo keyword in variable requires " + "thermo to use/init press"); + if (update->whichflag == 0) { + if (pressure->invoked_scalar != update->ntimestep) + error->all(FLERR,"Compute used in variable thermo keyword between runs " + "is not current"); + } else if (!(pressure->invoked_flag & INVOKED_SCALAR)) { + pressure->compute_scalar(); + pressure->invoked_flag |= INVOKED_SCALAR; + } + compute_enthalpy(); + + } else if (strcmp(word,"ecouple") == 0) compute_ecouple(); + } else if (strcmp(word,"econserve") == 0) { + if (!pe) + error->all(FLERR, + "Thermo keyword in variable requires thermo to use/init pe"); + if (update->whichflag == 0) { + if (pe->invoked_scalar != update->ntimestep) + error->all(FLERR,"Compute used in variable thermo keyword between runs " + "is not current"); + } else { + pe->compute_scalar(); + pe->invoked_flag |= INVOKED_SCALAR; + } + if (!temperature) + error->all(FLERR,"Thermo keyword in variable requires " + "thermo to use/init temp"); + if (update->whichflag == 0) { + if (temperature->invoked_scalar != update->ntimestep) + error->all(FLERR,"Compute used in variable thermo keyword between runs " + "is not current"); + } else if (!(temperature->invoked_flag & INVOKED_SCALAR)) { + temperature->compute_scalar(); + temperature->invoked_flag |= INVOKED_SCALAR; + } + compute_econserve(); + } else if (strcmp(word,"vol") == 0) compute_vol(); else if (strcmp(word,"density") == 0) compute_density(); else if (strcmp(word,"lx") == 0) compute_lx(); @@ -1721,27 +1756,28 @@ void Thermo::compute_ke() void Thermo::compute_etotal() { compute_pe(); - double ke = temperature->scalar; - ke *= 0.5 * temperature->dof * force->boltz; - if (normflag) ke /= natoms; - dvalue += ke; + double dvalue_pe = dvalue; + compute_ke(); + dvalue += dvalue_pe; } /* ---------------------------------------------------------------------- */ -void Thermo::compute_enthalpy() +void Thermo::compute_ecouple() { - compute_etotal(); - double etmp = dvalue; + dvalue = modify->ecouple(); +} - compute_vol(); - double vtmp = dvalue; - if (normflag) vtmp /= natoms; +/* ---------------------------------------------------------------------- */ - compute_press(); - double ptmp = dvalue; - - dvalue = etmp + ptmp*vtmp/(force->nktv2p); +void Thermo::compute_econserve() +{ + compute_pe(); + double dvalue_pe = dvalue; + compute_ke(); + double dvalue_ke = dvalue; + compute_ecouple(); + dvalue += dvalue_pe + dvalue_ke; } /* ---------------------------------------------------------------------- */ @@ -1867,6 +1903,24 @@ void Thermo::compute_etail() } else dvalue = 0.0; } +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_enthalpy() +{ + compute_etotal(); + double etmp = dvalue; + + compute_vol(); + double vtmp = dvalue; + if (normflag) vtmp /= natoms; + + compute_press(); + double ptmp = dvalue; + + dvalue = etmp + ptmp*vtmp/(force->nktv2p); +} + + /* ---------------------------------------------------------------------- */ void Thermo::compute_vol() diff --git a/src/thermo.h b/src/thermo.h index def8cead65..90ce41c6e4 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -143,7 +143,6 @@ class Thermo : protected Pointers { void compute_pe(); void compute_ke(); void compute_etotal(); - void compute_enthalpy(); void compute_evdwl(); void compute_ecoul(); @@ -156,6 +155,10 @@ class Thermo : protected Pointers { void compute_elong(); void compute_etail(); + void compute_enthalpy(); + void compute_ecouple(); + void compute_econserve(); + void compute_vol(); void compute_density(); void compute_lx();