diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt index f8e9f49073..23de1b1594 100644 --- a/doc/src/Section_commands.txt +++ b/doc/src/Section_commands.txt @@ -924,7 +924,7 @@ KOKKOS, o = USER-OMP, t = OPT. "tip4p/long (o)"_pair_coul.html, "tri/lj"_pair_tri_lj.html, "vashishta (o)"_pair_vashishta.html, -"vashishta/table (o)"_pair_vashishta_table.html, +"vashishta/table (o)"_pair_vashishta.html, "yukawa (go)"_pair_yukawa.html, "yukawa/colloid (go)"_pair_yukawa_colloid.html, "zbl (go)"_pair_zbl.html :tb(c=4,ea=c) diff --git a/doc/src/fix_dpd_energy.txt b/doc/src/fix_dpd_energy.txt index 9c4a55e964..db1fb0c703 100644 --- a/doc/src/fix_dpd_energy.txt +++ b/doc/src/fix_dpd_energy.txt @@ -27,9 +27,9 @@ the group at each timestep. It must be used in conjunction with a deterministic integrator (e.g. "fix nve"_fix_nve.html) that updates the particle positions and velocities. -For fix {dpd/energy}, the particle internal temperature is related to -the particle internal energy through a mesoparticle equation of state. -An additional fix must be specified that defines the equation of state +For fix {dpd/energy}, the particle internal temperature is related to +the particle internal energy through a mesoparticle equation of state. +An additional fix must be specified that defines the equation of state for each particle, e.g. "fix eos/cv"_fix_eos_cv.html. This fix must be used with the "pair_style @@ -37,7 +37,7 @@ dpd/fdt/energy"_pair_style.html command. Note that numerous variants of DPD can be specified by choosing an appropriate combination of the integrator and "pair_style -dpd/fdt/energy"_pair_style.html command. DPD under isoenergetic conditions +dpd/fdt/energy"_pair_style.html command. DPD under isoenergetic conditions can be specified by using fix {dpd/energy}, fix {nve} and pair_style {dpd/fdt/energy}. DPD under isoenthalpic conditions can be specified by using fix {dpd/energy}, fix {nph} and pair_style @@ -46,7 +46,7 @@ examples/USER/dpd directory. :line -[Restrictions:] +[Restrictions:] This command is part of the USER-DPD package. It is only enabled if LAMMPS was built with that package. See the "Making @@ -55,11 +55,11 @@ LAMMPS"_Section_start.html#start_3 section for more info. This fix must be used with an additional fix that specifies time integration, e.g. "fix nve"_fix_nve.html. -The fix {dpd/energy} requires the {dpd} "atom_style"_atom_style.html -to be used in order to properly account for the particle internal +The fix {dpd/energy} requires the {dpd} "atom_style"_atom_style.html +to be used in order to properly account for the particle internal energies and temperature. -The fix {dpd/energy} must be used with an additional fix that specifies the +The fix {dpd/energy} must be used with an additional fix that specifies the mesoparticle equation of state for each particle. [Related commands:] diff --git a/doc/src/fix_wall_gran.txt b/doc/src/fix_wall_gran.txt index f30872186d..40cd5c8c74 100644 --- a/doc/src/fix_wall_gran.txt +++ b/doc/src/fix_wall_gran.txt @@ -163,8 +163,8 @@ Any dimension (xyz) that has a granular wall must be non-periodic. [Related commands:] -"fix move"_fix_move.html, -"fix wall/gran/region"_fix_wall_gran_region.html, +"fix move"_fix_move.html, +"fix wall/gran/region"_fix_wall_gran_region.html, "pair_style granular"_pair_gran.html [Default:] none diff --git a/doc/src/fix_wall_gran_region.txt b/doc/src/fix_wall_gran_region.txt index 63092486eb..478f8b9b23 100644 --- a/doc/src/fix_wall_gran_region.txt +++ b/doc/src/fix_wall_gran_region.txt @@ -84,7 +84,7 @@ is up to you to ensure that the region location with respect to periodic or non-periodic boundaries is specified appropriately via the "region"_region.html and "boundary"_boundary.html commands when using a region as a wall that bounds particle motion. - + NOTE: For primitive regions with sharp corners and/or edges (e.g. a block or cylinder), wall/particle forces are computed accurately for both interior and exterior regions. For {union} and {intersect} @@ -190,7 +190,7 @@ LAMMPS"_Section_start.html#start_3 section for more info. [Related commands:] "fix_move"_fix_move.html, -"fix wall/gran"_fix_wall_gran.html, +"fix wall/gran"_fix_wall_gran.html, "fix wall/region"_fix_wall_region.html, "pair_style granular"_pair_gran.html, "region"_region.html diff --git a/doc/src/fixes.txt b/doc/src/fixes.txt index c8e989cc62..05a320085b 100644 --- a/doc/src/fixes.txt +++ b/doc/src/fixes.txt @@ -32,6 +32,7 @@ Fixes :h1 fix_drag fix_drude fix_drude_transform + fix_dpd_energy fix_dt_reset fix_efield fix_ehex @@ -148,6 +149,7 @@ Fixes :h1 fix_viscous fix_wall fix_wall_gran + fix_wall_gran_region fix_wall_piston fix_wall_reflect fix_wall_region diff --git a/doc/src/pair_dpd_fdt.txt b/doc/src/pair_dpd_fdt.txt index 86cd27b131..ffd138ddbd 100644 --- a/doc/src/pair_dpd_fdt.txt +++ b/doc/src/pair_dpd_fdt.txt @@ -33,12 +33,12 @@ pair_coeff * * 3.0 1.0 0.1 2.5 :pre [Description:] -Styles {dpd/fdt} and {dpd/fdt/energy} compute the force for dissipative -particle dynamics (DPD) simulations. The {dpd/fdt} style is used to -perform DPD simulations under isothermal and isobaric conditions, -while the {dpd/fdt/energy} style is used to perform DPD simulations -under isoenergetic and isoenthalpic conditions (see "(Lisal)"_#Lisal). -For DPD simulations in general, the force on atom I due to atom J is +Styles {dpd/fdt} and {dpd/fdt/energy} compute the force for dissipative +particle dynamics (DPD) simulations. The {dpd/fdt} style is used to +perform DPD simulations under isothermal and isobaric conditions, +while the {dpd/fdt/energy} style is used to perform DPD simulations +under isoenergetic and isoenthalpic conditions (see "(Lisal)"_#Lisal). +For DPD simulations in general, the force on atom I due to atom J is given as a sum of 3 terms :c,image(Eqs/pair_dpd.jpg) @@ -48,20 +48,20 @@ a random force. Rij is a unit vector in the direction Ri - Rj, Vij is the vector difference in velocities of the two atoms = Vi - Vj, alpha is a Gaussian random number with zero mean and unit variance, dt is the timestep size, and w(r) is a weighting factor that varies between -0 and 1. Rc is the cutoff. The weighting factor, omega_ij, varies +0 and 1. Rc is the cutoff. The weighting factor, omega_ij, varies between 0 and 1, and is chosen to have the following functional form: :c,image(Eqs/pair_dpd_omega.jpg) -Note that alternative definitions of the weighting function exist, but -would have to be implemented as a separate pair style command. +Note that alternative definitions of the weighting function exist, but +would have to be implemented as a separate pair style command. For style {dpd/fdt}, the fluctuation-dissipation theorem defines gamma to be set equal to sigma*sigma/(2 T), where T is the set point temperature specified as a pair style parameter in the above examples. -The following coefficients must be defined for each pair of atoms types -via the "pair_coeff"_pair_coeff.html command as in the examples above, -or in the data file or restart files read by the +The following coefficients must be defined for each pair of atoms types +via the "pair_coeff"_pair_coeff.html command as in the examples above, +or in the data file or restart files read by the "read_data"_read_data.html or "read_restart"_read_restart.html commands: A (force units) @@ -71,28 +71,28 @@ cutoff (distance units) :ul The last coefficient is optional. If not specified, the global DPD cutoff is used. -Style {dpd/fdt/energy} is used to perform DPD simulations -under isoenergetic and isoenthalpic conditions. The fluctuation-dissipation +Style {dpd/fdt/energy} is used to perform DPD simulations +under isoenergetic and isoenthalpic conditions. The fluctuation-dissipation theorem defines gamma to be set equal to sigma*sigma/(2 dpdTheta), where -dpdTheta is the average internal temperature for the pair. The particle -internal temperature is related to the particle internal energy through -a mesoparticle equation of state (see "fix eos"_fix.html). The -differential internal conductive and mechanical energies are computed +dpdTheta is the average internal temperature for the pair. The particle +internal temperature is related to the particle internal energy through +a mesoparticle equation of state (see "fix eos"_fix.html). The +differential internal conductive and mechanical energies are computed within style {dpd/fdt/energy} as: :c,image(Eqs/pair_dpd_energy.jpg) -where +where :c,image(Eqs/pair_dpd_energy_terms.jpg) -Zeta_ij^q is a second Gaussian random number with zero mean and unit -variance that is used to compute the internal conductive energy. The +Zeta_ij^q is a second Gaussian random number with zero mean and unit +variance that is used to compute the internal conductive energy. The fluctuation-dissipation theorem defines alpha*alpha to be set equal to 2*kB*kappa, where kappa is the mesoparticle thermal -conductivity parameter. The following coefficients must be defined for -each pair of atoms types via the "pair_coeff"_pair_coeff.html -command as in the examples above, or in the data file or restart files +conductivity parameter. The following coefficients must be defined for +each pair of atoms types via the "pair_coeff"_pair_coeff.html +command as in the examples above, or in the data file or restart files read by the "read_data"_read_data.html or "read_restart"_read_restart.html commands: @@ -109,18 +109,18 @@ The pairwise energy associated with styles {dpd/fdt} and shifted to be zero at the cutoff distance Rc. The pairwise virial is calculated using only the conservative term. -The forces computed through the {dpd/fdt} and {dpd/fdt/energy} styles +The forces computed through the {dpd/fdt} and {dpd/fdt/energy} styles can be integrated with the velocity-Verlet integration scheme or the -Shardlow splitting integration scheme described by "(Lisal)"_#Lisal. -In the cases when these pair styles are combined with the -"fix shardlow"_fix_shardlow.html, these pair styles differ from the +Shardlow splitting integration scheme described by "(Lisal)"_#Lisal. +In the cases when these pair styles are combined with the +"fix shardlow"_fix_shardlow.html, these pair styles differ from the other dpd styles in that the dissipative and random forces are split -from the force calculation and are not computed within the pair style. -Thus, only the conservative force is computed by the pair style, -while the stochastic integration of the dissipative and random forces -are handled through the Shardlow splitting algorithm approach. The -Shardlow splitting algorithm is advantageous, especially when -performing DPD under isoenergetic conditions, as it allows +from the force calculation and are not computed within the pair style. +Thus, only the conservative force is computed by the pair style, +while the stochastic integration of the dissipative and random forces +are handled through the Shardlow splitting algorithm approach. The +Shardlow splitting algorithm is advantageous, especially when +performing DPD under isoenergetic conditions, as it allows significantly larger timesteps to be taken. :line @@ -132,7 +132,7 @@ enabled if LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#start_3 section for more info. Pair styles {dpd/fdt} and {dpd/fdt/energy} require use of the -"communicate vel yes"_communicate.html option so that velocites are +"comm_modify vel yes"_comm_modify.html option so that velocites are stored by ghost atoms. Pair style {dpd/fdt/energy} requires "atom_style dpd"_atom_style.html diff --git a/src/.gitignore b/src/.gitignore index 857ac769b4..164ad0f119 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -489,6 +489,8 @@ /fix_wall_colloid.h /fix_wall_gran.cpp /fix_wall_gran.h +/fix_wall_gran_region.cpp +/fix_wall_gran_region.h /fix_wall_piston.cpp /fix_wall_piston.h /fix_wall_srd.cpp diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index 3ed7cf3d90..96b1f5071f 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -15,9 +15,9 @@ Contributing authors: Dan Bolintineanu (SNL) ------------------------------------------------------------------------- */ -#include "math.h" -#include "stdlib.h" -#include "string.h" +#include +#include +#include #include "fix_wall_gran_region.h" #include "region.h" #include "atom.h" diff --git a/src/USER-DPD/fix_dpd_energy.cpp b/src/USER-DPD/fix_dpd_energy.cpp index e6b7a9e83e..05907a5fcf 100644 --- a/src/USER-DPD/fix_dpd_energy.cpp +++ b/src/USER-DPD/fix_dpd_energy.cpp @@ -11,8 +11,8 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "stdio.h" -#include "string.h" +#include +#include #include "fix_dpd_energy.h" #include "atom.h" #include "force.h" diff --git a/src/USER-DPD/fix_eos_table_rx.cpp b/src/USER-DPD/fix_eos_table_rx.cpp index 4cd6ffc755..85c3a6dd2a 100644 --- a/src/USER-DPD/fix_eos_table_rx.cpp +++ b/src/USER-DPD/fix_eos_table_rx.cpp @@ -687,7 +687,7 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai) double maxit = 100; double temp; double delta = 0.001; - + // Store the current thetai in t1 t1 = MAX(thetai,tb->lo); t1 = MIN(t1,tb->hi); @@ -731,7 +731,7 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai) if(it==maxit){ char str[256]; sprintf(str,"Maxit exceeded in secant solver: id=%d ui=%lf thetai=%lf t1=%lf t2=%lf f1=%lf f2=%lf\n",id,ui,thetai,t1,t2,f1,f2); - if(isnan(f1) || isnan(f2) || isnan(ui) || isnan(thetai) || isnan(t1) || isnan(t2)) + if(isnan(f1) || isnan(f2) || isnan(ui) || isnan(thetai) || isnan(t1) || isnan(t2)) error->one(FLERR,"NaN detected in secant solver."); error->one(FLERR,str); } diff --git a/src/USER-DPD/fix_rx.cpp b/src/USER-DPD/fix_rx.cpp index 96f7f44422..67b134e53e 100644 --- a/src/USER-DPD/fix_rx.cpp +++ b/src/USER-DPD/fix_rx.cpp @@ -231,7 +231,7 @@ FixRX::~FixRX() memory->destroy( sparseKinetics_inu ); memory->destroy( sparseKinetics_isIntegralReaction ); } -} +} /* ---------------------------------------------------------------------- */ @@ -1720,7 +1720,7 @@ void FixRX::computeLocalTemperature() // Lucy's Weight Function if(wtFlag==LUCY){ - wij = (1.0+3.0*ratio) * (1.0-ratio)*(1.0-ratio)*(1.0-ratio); + wij = (1.0+3.0*ratio) * (1.0-ratio)*(1.0-ratio)*(1.0-ratio); dpdThetaLocal[i] += wij/dpdTheta[j]; if (newton_pair || j < nlocal) dpdThetaLocal[j] += wij/dpdTheta[i]; diff --git a/src/USER-DPD/pair_exp6_rx.cpp b/src/USER-DPD/pair_exp6_rx.cpp index 060ef2d7c1..ad2dbf1465 100644 --- a/src/USER-DPD/pair_exp6_rx.cpp +++ b/src/USER-DPD/pair_exp6_rx.cpp @@ -451,7 +451,7 @@ void PairExp6rx::compute(int eflag, int vflag) // // Apply Mixing Rule to get the overall force for the CG pair // - if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairOldEXP6_12; + if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairOldEXP6_12; else fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairOldEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*fpairOldEXP6_21; f[i][0] += delx*fpair; @@ -588,7 +588,7 @@ void PairExp6rx::coeff(int narg, char **arg) { // Set isite1 and isite2 parameters based on site1 and site2 strings. - + if (strcmp(site1,"1fluid") == 0) isite1 = oneFluidApproxParameter; else @@ -602,7 +602,7 @@ void PairExp6rx::coeff(int narg, char **arg) else isite1 = isp; } - + if (strcmp(site2,"1fluid") == 0) isite2 = oneFluidApproxParameter; else @@ -616,7 +616,7 @@ void PairExp6rx::coeff(int narg, char **arg) else isite2 = isp; } - + // Set the interaction potential type to the enumerated type. for (int iparam = 0; iparam < nparams; ++iparam) { diff --git a/src/USER-DPD/pair_multi_lucy_rx.cpp b/src/USER-DPD/pair_multi_lucy_rx.cpp index de86733725..cf10a70756 100644 --- a/src/USER-DPD/pair_multi_lucy_rx.cpp +++ b/src/USER-DPD/pair_multi_lucy_rx.cpp @@ -234,7 +234,7 @@ void PairMultiLucyRX::compute(int eflag, int vflag) } else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx"); - if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair; + if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair; else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair; fx_i += delx*fpair; @@ -935,7 +935,7 @@ void PairMultiLucyRX::getParams(int id, double &fractionOld1, double &fractionOl nTotal = 0.0; nTotalOld = 0.0; for (int ispecies = 0; ispecies < nspecies; ispecies++){ - nTotal += atom->dvector[ispecies][id]; + nTotal += atom->dvector[ispecies][id]; nTotalOld += atom->dvector[ispecies+nspecies][id]; } diff --git a/src/USER-DPD/pair_table_rx.cpp b/src/USER-DPD/pair_table_rx.cpp index 44a9d76022..74b9222a0b 100644 --- a/src/USER-DPD/pair_table_rx.cpp +++ b/src/USER-DPD/pair_table_rx.cpp @@ -186,7 +186,7 @@ void PairTableRX::compute(int eflag, int vflag) value = tb->f[itable] + fraction*tb->df[itable]; fpair = factor_lj * value; } - if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair; + if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair; else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair; fx_i += delx*fpair; @@ -1102,7 +1102,7 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq, fforce = factor_lj * value; } - if (isite1 == isite2) fforce = sqrt(fraction1_i*fraction2_j)*fforce; + if (isite1 == isite2) fforce = sqrt(fraction1_i*fraction2_j)*fforce; else fforce = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*fforce; if (tabstyle == LOOKUP) @@ -1146,17 +1146,17 @@ void PairTableRX::getParams(int id, double &fractionOld1, double &fractionOld2, double nTotal = 0.0; double nTotalOld = 0.0; for (int ispecies = 0; ispecies < nspecies; ++ispecies){ - nTotal += atom->dvector[ispecies][id]; + nTotal += atom->dvector[ispecies][id]; nTotalOld += atom->dvector[ispecies+nspecies][id]; } if(nTotal < 1e-8 || nTotalOld < 1e-8) error->all(FLERR,"The number of molecules in CG particle is less than 1e-8."); - if (isOneFluid(isite1) == false){ + if (isOneFluid(isite1) == false){ fractionOld1 = atom->dvector[isite1+nspecies][id]/nTotalOld; fraction1 = atom->dvector[isite1][id]/nTotal; } - if (isOneFluid(isite2) == false){ + if (isOneFluid(isite2) == false){ fractionOld2 = atom->dvector[isite2+nspecies][id]/nTotalOld; fraction2 = atom->dvector[isite2][id]/nTotal; } diff --git a/src/timer.cpp b/src/timer.cpp index 8711231889..583c6a2b3f 100644 --- a/src/timer.cpp +++ b/src/timer.cpp @@ -98,6 +98,7 @@ Timer::Timer(LAMMPS *lmp) : Pointers(lmp) _timeout = -1.0; _checkfreq = 10; _nextcheck = -1; + _laststep = -1; this->_stamp(RESET); } @@ -215,6 +216,7 @@ void Timer::set_wall(enum ttype which, double newtime) void Timer::init_timeout() { + _laststep = -1; if (_timeout < 0) _nextcheck = -1; else @@ -247,6 +249,14 @@ void Timer::print_timeout(FILE *fp) /* ---------------------------------------------------------------------- */ +void Timer::force_timeout() +{ + _timeout = 0.0; + _nextcheck = _laststep + 1; +} + +/* ---------------------------------------------------------------------- */ + bool Timer::_check_timeout() { double walltime = MPI_Wtime() - timeout_start; diff --git a/src/timer.h b/src/timer.h index 903a182128..ff9c66db7c 100644 --- a/src/timer.h +++ b/src/timer.h @@ -64,7 +64,7 @@ class Timer : protected Pointers { void init_timeout(); // trigger enforced timeout - void force_timeout() { _timeout = 0.0; }; + void force_timeout(); // get remaining time in seconds. 0.0 if inactive, negative if expired double get_timeout_remain(); @@ -75,6 +75,7 @@ class Timer : protected Pointers { // check for timeout. inline wrapper around internal // function to reduce overhead in case there is no check. bool check_timeout(int step) { + _laststep = step; if (_nextcheck != step) return false; else return _check_timeout(); } @@ -91,7 +92,8 @@ class Timer : protected Pointers { int _sync; // if nonzero, synchronize tasks before setting the timer int _timeout; // max allowed wall time in seconds. infinity if negative int _checkfreq; // frequency of timeout checking - int _nextcheck; // timestep number of next timeout check + int _nextcheck; // loop number of next timeout check + int _laststep; // loop number of last iteration // update one specific timer array void _stamp(enum ttype);