ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #217 from akohlmey/small-fixes

Browse Source 
Collected small changes and bugfixes
This commit is contained in:
sjplimp
2016-10-12 07:32:08 -06:00
committed by GitHub
parent 4e3a55047f f8a26dd158
commit 4a5d9eaae2
16 changed files with 86 additions and 70 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Section_commands.txt
View File

@ -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)

16
doc/src/fix_dpd_energy.txt
View File

@ -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:]

4
doc/src/fix_wall_gran.txt
View File

@ -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

4
doc/src/fix_wall_gran_region.txt
View File

@ -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

2
doc/src/fixes.txt
View File

@ -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

70
doc/src/pair_dpd_fdt.txt
View File

@ -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

2
src/.gitignore vendored
View File

@ -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

6
src/GRANULAR/fix_wall_gran_region.cpp
View File

@ -15,9 +15,9 @@
Contributing authors: Dan Bolintineanu (SNL)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "fix_wall_gran_region.h"
#include "region.h"
#include "atom.h"

4
src/USER-DPD/fix_dpd_energy.cpp
View File

@ -11,8 +11,8 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "stdio.h"
#include "string.h"
#include <stdio.h>
#include <string.h>
#include "fix_dpd_energy.h"
#include "atom.h"
#include "force.h"

4
src/USER-DPD/fix_eos_table_rx.cpp
View File

@ -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);
}

4
src/USER-DPD/fix_rx.cpp
View File

@ -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];

8
src/USER-DPD/pair_exp6_rx.cpp
View File

@ -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)
{

4
src/USER-DPD/pair_multi_lucy_rx.cpp
View File

@ -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];
}

10
src/USER-DPD/pair_table_rx.cpp
View File

@ -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;
}

10
src/timer.cpp
View File

@ -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;

6
src/timer.h
View File

@ -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);