diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index b3c0a70e4a..a334ba0368 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -171,8 +171,8 @@ set(DEFAULT_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS DIPOLE GRANULAR USER-BOCS USER-CGDNA USER-MESO USER-CGSDK USER-COLVARS USER-DIFFRACTION USER-DPD USER-DRUDE USER-EFF USER-FEP USER-H5MD USER-LB USER-MANIFOLD USER-MEAMC USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF - USER-PHONON USER-PTM USER-QTB USER-REAXC USER-SCAFACOS USER-SMD USER-SMTBQ - USER-SPH USER-TALLY USER-UEF USER-VTK USER-QUIP USER-QMMM) + USER-PHONON USER-PTM USER-QTB USER-REAXC USER-SCAFACOS USER-SDPD USER-SMD + USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK USER-QUIP USER-QMMM) set(ACCEL_PACKAGES USER-OMP KOKKOS OPT USER-INTEL GPU) set(OTHER_PACKAGES CORESHELL QEQ) foreach(PKG ${DEFAULT_PACKAGES}) diff --git a/cmake/presets/all_off.cmake b/cmake/presets/all_off.cmake index f7e90ddbb4..9021dedc08 100644 --- a/cmake/presets/all_off.cmake +++ b/cmake/presets/all_off.cmake @@ -8,7 +8,7 @@ set(USER_PACKAGES USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-CGSDK USER-COLVA USER-INTEL USER-LB USER-MANIFOLD USER-MEAMC USER-MESO USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF USER-OMP USER-PHONON USER-QMMM USER-QTB - USER-QUIP USER-REAXC USER-SMD USER-SMTBQ USER-SPH USER-TALLY + USER-QUIP USER-REAXC USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK) set(PACKAGES_WITH_LIB COMPRESS GPU KIM KOKKOS LATTE MEAM MPIIO MSCG POEMS PYTHON REAX VORONOI diff --git a/cmake/presets/all_on.cmake b/cmake/presets/all_on.cmake index 2c6f67904e..ef17e38914 100644 --- a/cmake/presets/all_on.cmake +++ b/cmake/presets/all_on.cmake @@ -8,7 +8,7 @@ set(USER_PACKAGES USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-CGSDK USER-COLVA USER-INTEL USER-LB USER-MANIFOLD USER-MEAMC USER-MESO USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF USER-OMP USER-PHONON USER-QMMM USER-QTB - USER-QUIP USER-REAXC USER-SMD USER-SMTBQ USER-SPH USER-TALLY + USER-QUIP USER-REAXC USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK) set(PACKAGES_WITH_LIB COMPRESS GPU KIM KOKKOS LATTE MEAM MPIIO MSCG POEMS PYTHON REAX VORONOI diff --git a/cmake/presets/manual_selection.cmake b/cmake/presets/manual_selection.cmake index 6c03d983f6..35a45a5c87 100644 --- a/cmake/presets/manual_selection.cmake +++ b/cmake/presets/manual_selection.cmake @@ -61,6 +61,7 @@ set(PKG_USER-QMMM OFF CACHE BOOL "" FORCE) set(PKG_USER-QTB OFF CACHE BOOL "" FORCE) set(PKG_USER-QUIP OFF CACHE BOOL "" FORCE) set(PKG_USER-REAXC OFF CACHE BOOL "" FORCE) +set(PKG_USER-SDPD OFF CACHE BOOL "" FORCE) set(PKG_USER-SMD OFF CACHE BOOL "" FORCE) set(PKG_USER-SMTBQ OFF CACHE BOOL "" FORCE) set(PKG_USER-SPH OFF CACHE BOOL "" FORCE) diff --git a/cmake/presets/nolib.cmake b/cmake/presets/nolib.cmake index cd603aa804..54db12a851 100644 --- a/cmake/presets/nolib.cmake +++ b/cmake/presets/nolib.cmake @@ -8,7 +8,7 @@ set(USER_PACKAGES USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-CGSDK USER-COLVA USER-INTEL USER-LB USER-MANIFOLD USER-MEAMC USER-MESO USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF USER-OMP USER-PHONON USER-QMMM USER-QTB - USER-QUIP USER-REAXC USER-SMD USER-SMTBQ USER-SPH USER-TALLY + USER-QUIP USER-REAXC USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK) set(PACKAGES_WITH_LIB COMPRESS GPU KIM KOKKOS LATTE MEAM MPIIO MSCG POEMS PYTHON REAX VORONOI diff --git a/cmake/presets/std.cmake b/cmake/presets/std.cmake index 36da897957..4176aba44e 100644 --- a/cmake/presets/std.cmake +++ b/cmake/presets/std.cmake @@ -8,7 +8,7 @@ set(USER_PACKAGES USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-CGSDK USER-COLVA USER-INTEL USER-LB USER-MANIFOLD USER-MEAMC USER-MESO USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF USER-OMP USER-PHONON USER-QMMM USER-QTB - USER-QUIP USER-REAXC USER-SMD USER-SMTBQ USER-SPH USER-TALLY + USER-QUIP USER-REAXC USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK) set(PACKAGES_WITH_LIB COMPRESS GPU KIM KOKKOS LATTE MEAM MPIIO MSCG POEMS PYTHON REAX VORONOI diff --git a/cmake/presets/std_nolib.cmake b/cmake/presets/std_nolib.cmake index 9bffefcbe0..aa067f2ba0 100644 --- a/cmake/presets/std_nolib.cmake +++ b/cmake/presets/std_nolib.cmake @@ -8,7 +8,7 @@ set(USER_PACKAGES USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-CGSDK USER-COLVA USER-INTEL USER-LB USER-MANIFOLD USER-MEAMC USER-MESO USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF USER-OMP USER-PHONON USER-QMMM USER-QTB - USER-QUIP USER-REAXC USER-SMD USER-SMTBQ USER-SPH USER-TALLY + USER-QUIP USER-REAXC USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK) set(PACKAGES_WITH_LIB COMPRESS GPU KIM KOKKOS LATTE MEAM MPIIO MSCG POEMS PYTHON REAX VORONOI diff --git a/cmake/presets/user.cmake b/cmake/presets/user.cmake index cb81b67558..afb030d3d6 100644 --- a/cmake/presets/user.cmake +++ b/cmake/presets/user.cmake @@ -8,7 +8,7 @@ set(USER_PACKAGES USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-CGSDK USER-COLVA USER-INTEL USER-LB USER-MANIFOLD USER-MEAMC USER-MESO USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF USER-OMP USER-PHONON USER-QMMM USER-QTB - USER-QUIP USER-REAXC USER-SMD USER-SMTBQ USER-SPH USER-TALLY + USER-QUIP USER-REAXC USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK) set(PACKAGES_WITH_LIB COMPRESS GPU KIM KOKKOS LATTE MEAM MPIIO MSCG POEMS PYTHON REAX VORONOI diff --git a/doc/src/Commands_fix.txt b/doc/src/Commands_fix.txt index 4144c7eb90..77684390bd 100644 --- a/doc/src/Commands_fix.txt +++ b/doc/src/Commands_fix.txt @@ -94,6 +94,7 @@ OPT. "lineforce"_fix_lineforce.html, "manifoldforce"_fix_manifoldforce.html, "meso"_fix_meso.html, +"meso/move"_fix_meso_move.html, "meso/stationary"_fix_meso_stationary.html, "momentum (k)"_fix_momentum.html, "move"_fix_move.html, @@ -172,6 +173,7 @@ OPT. "restrain"_fix_restrain.html, "rhok"_fix_rhok.html, "rigid (o)"_fix_rigid.html, +"rigid/meso"_fix_rigid_meso.html, "rigid/nph (o)"_fix_rigid.html, "rigid/nph/small"_fix_rigid.html, "rigid/npt (o)"_fix_rigid.html, diff --git a/doc/src/Commands_pair.txt b/doc/src/Commands_pair.txt index a74574367c..7af2cc9bae 100644 --- a/doc/src/Commands_pair.txt +++ b/doc/src/Commands_pair.txt @@ -198,6 +198,7 @@ OPT. "reax/c (ko)"_pair_reaxc.html, "rebo (io)"_pair_airebo.html, "resquared (go)"_pair_resquared.html, +"sdpd/taitwater/isothermal"_pair_sdpd_taitwater_isothermal.html, "smd/hertz"_pair_smd_hertz.html, "smd/tlsph"_pair_smd_tlsph.html, "smd/tri_surface"_pair_smd_triangulated_surface.html, diff --git a/doc/src/Packages_details.txt b/doc/src/Packages_details.txt index 031858e846..be9006316a 100644 --- a/doc/src/Packages_details.txt +++ b/doc/src/Packages_details.txt @@ -95,6 +95,7 @@ as contained in the file name. "USER-QUIP"_#PKG-USER-QUIP, "USER-REAXC"_#PKG-USER-REAXC, "USER-SCAFACOS"_#PKG-USER-SCAFACOS, +"USER-SDPD"_#PKG-USER-SDPD, "USER-SMD"_#PKG-USER-SMD, "USER-SMTBQ"_#PKG-USER-SMTBQ, "USER-SPH"_#PKG-USER-SPH, @@ -1916,6 +1917,31 @@ examples/USER/scafacos :ul :line +USER-SDPD package :link(PKG-USER-SDPD),h4 + +[Contents:] + +A pair style for smoothed dissipative particle dynamics (SDPD), which +is an extension of smoothed particle hydrodynamics (SPH) to mesoscale +where thermal fluctuations are important (see the +"USER-SPH package"_#PKG-USER-SPH). +Also two fixes for moving and rigid body integration of SPH/SDPD particles +(particles of atom_style meso). + +[Author:] Morteza Jalalvand (Institute for Advanced Studies in Basic +Sciences, Iran). + +[Supporting info:] + +src/USER-SDPD: filenames -> commands +src/USER-SDPD/README +"pair_style sdpd/taitwater/isothermal"_pair_sdpd_taitwater_isothermal.html +"fix meso/move"_fix_meso_move.html +"fix rigid/meso"_fix_rigid_meso.html +examples/USER/sdpd :ul + +:line + USER-SMD package :link(PKG-USER-SMD),h4 [Contents:] diff --git a/doc/src/Packages_user.txt b/doc/src/Packages_user.txt index 9bb3fbd18c..e1269b0d0c 100644 --- a/doc/src/Packages_user.txt +++ b/doc/src/Packages_user.txt @@ -68,6 +68,7 @@ Package, Description, Doc page, Example, Library "USER-QUIP"_Packages_details.html#PKG-USER-QUIP, QUIP/libatoms interface,"pair_style quip"_pair_quip.html, USER/quip, ext "USER-REAXC"_Packages_details.html#PKG-USER-REAXC, ReaxFF potential (C/C++) ,"pair_style reaxc"_pair_reaxc.html, reax, no "USER-SCAFACOS"_Packages_details.html#PKG-USER-SCAFACOS, wrapper on ScaFaCoS solver,"kspace_style scafacos"_kspace_style.html, USER/scafacos, ext +"USER-SDPD"_Packages_details.html#PKG-USER-SDPD, smoothed dissipative particle dynamics,"pair_style sdpd/taitwater/isothermal"_pair_sdpd_taitwater_isothermal, USER/sdpd, no "USER-SMD"_Packages_details.html#PKG-USER-SMD, smoothed Mach dynamics,"SMD User Guide"_PDF/SMD_LAMMPS_userguide.pdf, USER/smd, ext "USER-SMTBQ"_Packages_details.html#PKG-USER-SMTBQ, second moment tight binding QEq potential,"pair_style smtbq"_pair_smtbq.html, USER/smtbq, no "USER-SPH"_Packages_details.html#PKG-USER-SPH, smoothed particle hydrodynamics,"SPH User Guide"_PDF/SPH_LAMMPS_userguide.pdf, USER/sph, no diff --git a/doc/src/fix.txt b/doc/src/fix.txt index 27cc0467e2..8d00fa987d 100644 --- a/doc/src/fix.txt +++ b/doc/src/fix.txt @@ -237,6 +237,7 @@ accelerated styles exist. "lineforce"_fix_lineforce.html - constrain atoms to move in a line "manifoldforce"_fix_manifoldforce.html - "meso"_fix_meso.html - +"meso"_fix_meso_move.html - move mesoscopic SPH/SDPD particles in a prescribed fashion "meso/stationary"_fix_meso_stationary.html - "momentum"_fix_momentum.html - zero the linear and/or angular momentum of a group of atoms "move"_fix_move.html - move atoms in a prescribed fashion @@ -331,6 +332,7 @@ accelerated styles exist. "rigid/small/npt"_fix_rigid.html - constrain many small clusters of atoms to move as a rigid body with NPT integration "rigid/small/nve"_fix_rigid.html - constrain many small clusters of atoms to move as a rigid body with alternate NVE integration "rigid/small/nvt"_fix_rigid.html - constrain many small clusters of atoms to move as a rigid body with NVT integration +"rigid/meso"_fix_rigid_meso.html - constrain clusters of mesoscopic SPH/SDPD particles to move as a rigid body "rx"_fix_rx.html - "saed/vtk"_fix_saed_vtk.html - "setforce"_fix_setforce.html - set the force on each atom diff --git a/doc/src/fix_meso_move.txt b/doc/src/fix_meso_move.txt new file mode 100644 index 0000000000..0a222e79ed --- /dev/null +++ b/doc/src/fix_meso_move.txt @@ -0,0 +1,233 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +fix meso/move command :h3 + +[Syntax:] + +fix ID group-ID meso/move style args keyword values ... :pre + +ID, group-ID are documented in "fix"_fix.html command :ulb,l +meso/move = style name of this fix command :l +style = {linear} or {wiggle} or {rotate} or {variable} :l + {linear} args = Vx Vy Vz + Vx,Vy,Vz = components of velocity vector (velocity units), any component can be specified as NULL + {wiggle} args = Ax Ay Az period + Ax,Ay,Az = components of amplitude vector (distance units), any component can be specified as NULL + period = period of oscillation (time units) + {rotate} args = Px Py Pz Rx Ry Rz period + Px,Py,Pz = origin point of axis of rotation (distance units) + Rx,Ry,Rz = axis of rotation vector + period = period of rotation (time units) + {variable} args = v_dx v_dy v_dz v_vx v_vy v_vz + v_dx,v_dy,v_dz = 3 variable names that calculate x,y,z displacement as function of time, any component can be specified as NULL + v_vx,v_vy,v_vz = 3 variable names that calculate x,y,z velocity as function of time, any component can be specified as NULL :pre + +zero or more keyword/value pairs may be appended :l +keyword = {units} :l + {units} value = {box} or {lattice} :pre +:ule + +[Examples:] + +fix 1 boundary meso/move wiggle 3.0 0.0 0.0 1.0 units box +fix 2 boundary meso/move rotate 0.0 0.0 0.0 0.0 0.0 1.0 5.0 +fix 2 boundary meso/move variable v_myx v_myy NULL v_VX v_VY NULL :pre + +[Description:] + +Perform updates of position, velocity, internal energy and local +density for mesoscopic particles in the group each timestep using the +specified settings or formulas, without regard to forces on the +particles. This can be useful for boundary, solid bodies or other +particles, whose movement can influence nearby particles. + +The operation of this fix is exactly like that described by the +"fix move"_fix_move.html command, except that particles' density, +internal energy and extrapolated velocity are also updated. + +NOTE: The particles affected by this fix should not be time integrated +by other fixes (e.g. "fix meso"_fix_meso.html, "fix +meso/stationary"_fix_meso_stationary.html), since that will change their +positions and velocities twice. + +NOTE: As particles move due to this fix, they will pass thru periodic +boundaries and be remapped to the other side of the simulation box, +just as they would during normal time integration (e.g. via the "fix +meso"_fix_meso.html command). It is up to you to decide whether periodic +boundaries are appropriate with the kind of particle motion you are +prescribing with this fix. + +NOTE: As dicsussed below, particles are moved relative to their initial +position at the time the fix is specified. These initial coordinates +are stored by the fix in "unwrapped" form, by using the image flags +associated with each particle. See the "dump custom"_dump.html command +for a discussion of "unwrapped" coordinates. See the Atoms section of +the "read_data"_read_data.html command for a discussion of image flags +and how they are set for each particle. You can reset the image flags +(e.g. to 0) before invoking this fix by using the "set image"_set.html +command. + +:line + +The {linear} style moves particles at a constant velocity, so that their +position {X} = (x,y,z) as a function of time is given in vector +notation as + +X(t) = X0 + V * delta :pre + +where {X0} = (x0,y0,z0) is their position at the time the fix is +specified, {V} is the specified velocity vector with components +(Vx,Vy,Vz), and {delta} is the time elapsed since the fix was +specified. This style also sets the velocity of each particle to V = +(Vx,Vy,Vz). If any of the velocity components is specified as NULL, +then the position and velocity of that component is time integrated +the same as the "fix meso"_fix_meso.html command would perform, using +the corresponding force component on the particle. + +Note that the {linear} style is identical to using the {variable} +style with an "equal-style variable"_variable.html that uses the +vdisplace() function. E.g. + +variable V equal 10.0 +variable x equal vdisplace(0.0,$V) +fix 1 boundary move variable v_x NULL NULL v_V NULL NULL :pre + +The {wiggle} style moves particles in an oscillatory fashion, so that +their position {X} = (x,y,z) as a function of time is given in vector +notation as + +X(t) = X0 + A sin(omega*delta) :pre + +where {X0} = (x0,y0,z0) is their position at the time the fix is +specified, {A} is the specified amplitude vector with components +(Ax,Ay,Az), {omega} is 2 PI / {period}, and {delta} is the time +elapsed since the fix was specified. This style also sets the +velocity of each particle to the time derivative of this expression. +If any of the amplitude components is specified as NULL, then the +position and velocity of that component is time integrated the same as +the "fix meso"_fix_meso.html command would perform, using the +corresponding force component on the particle. + +Note that the {wiggle} style is identical to using the {variable} +style with "equal-style variables"_variable.html that use the +swiggle() and cwiggle() functions. E.g. + +variable A equal 10.0 +variable T equal 5.0 +variable omega equal 2.0*PI/$T +variable x equal swiggle(0.0,$A,$T) +variable v equal v_omega*($A-cwiggle(0.0,$A,$T)) +fix 1 boundary move variable v_x NULL NULL v_v NULL NULL :pre + +The {rotate} style rotates particles around a rotation axis {R} = +(Rx,Ry,Rz) that goes thru a point {P} = (Px,Py,Pz). The {period} of +the rotation is also specified. The direction of rotation for the +particles around the rotation axis is consistent with the right-hand +rule: if your right-hand thumb points along {R}, then your fingers wrap +around the axis in the direction of rotation. + +This style also sets the velocity of each particle to (omega cross +Rperp) where omega is its angular velocity around the rotation axis and +Rperp is a perpendicular vector from the rotation axis to the particle. + +The {variable} style allows the position and velocity components of +each particle to be set by formulas specified via the +"variable"_variable.html command. Each of the 6 variables is +specified as an argument to the fix as v_name, where name is the +variable name that is defined elsewhere in the input script. + +Each variable must be of either the {equal} or {atom} style. +{Equal}-style variables compute a single numeric quantity, that can be +a function of the timestep as well as of other simulation values. +{Atom}-style variables compute a numeric quantity for each particle, that +can be a function per-atom quantities, such as the particle's position, as +well as of the timestep and other simulation values. Note that this +fix stores the original coordinates of each particle (see note below) so +that per-atom quantity can be used in an atom-style variable formula. +See the "variable"_variable.html command for details. + +The first 3 variables (v_dx,v_dy,v_dz) specified for the {variable} +style are used to calculate a displacement from the particle's original +position at the time the fix was specified. The second 3 variables +(v_vx,v_vy,v_vz) specified are used to compute a velocity for each +particle. + +Any of the 6 variables can be specified as NULL. If both the +displacement and velocity variables for a particular x,y,z component +are specified as NULL, then the position and velocity of that +component is time integrated the same as the "fix meso"_fix_meso.html +command would perform, using the corresponding force component on the +particle. If only the velocity variable for a component is specified as +NULL, then the displacement variable will be used to set the position +of the particle, and its velocity component will not be changed. If only +the displacement variable for a component is specified as NULL, then +the velocity variable will be used to set the velocity of the particle, +and the position of the particle will be time integrated using that +velocity. + +The {units} keyword determines the meaning of the distance units used +to define the {linear} velocity and {wiggle} amplitude and {rotate} +origin. This setting is ignored for the {variable} style. A {box} +value selects standard units as defined by the "units"_units.html +command, e.g. velocity in Angstroms/fmsec and amplitude and position +in Angstroms for units = real. A {lattice} value means the velocity +units are in lattice spacings per time and the amplitude and position +are in lattice spacings. The "lattice"_lattice.html command must have +been previously used to define the lattice spacing. Each of these 3 +quantities may be dependent on the x,y,z dimension, since the lattice +spacings can be different in x,y,z. + +:line + +[Restart, fix_modify, output, run start/stop, minimize info:] + +This fix writes the original coordinates of moving particles to "binary +restart files"_restart.html, as well as the initial timestep, so that +the motion can be continuous in a restarted simulation. See the +"read_restart"_read_restart.html command for info on how to re-specify +a fix in an input script that reads a restart file, so that the +operation of the fix continues in an uninterrupted fashion. + +NOTE: Because the move positions are a function of the current +timestep and the initial timestep, you cannot reset the timestep to a +different value after reading a restart file, if you expect a fix move +command to work in an uninterrupted fashion. + +None of the "fix_modify"_fix_modify.html options are relevant to this +fix. + +This fix produces a per-atom array which can be accessed by various +"output commands"_Howto_output.html. The number of columns for each +atom is 3, and the columns store the original unwrapped x,y,z coords +of each particle. The per-atom values can be accessed on any timestep. + +No parameter of this fix can be used with the {start/stop} keywords of +the "run"_run.html command. + +This fix is not invoked during "energy minimization"_minimize.html. + +[Restrictions:] + +This fix is part of the USER-SDPD package. It is only enabled if +LAMMPS was built with that package. See the "Build +package"_Build_package.html doc page for more info. + +This fix requires that atoms store density and internal energy as +defined by the "atom_style meso"_atom_style.html command. + +All particles in the group must be mesoscopic SPH/SDPD particles. + +[Related commands:] + +"fix move"_fix_move.html, "fix meso"_fix_meso.html, +"displace_atoms"_displace_atoms.html + +[Default:] + +The option default is units = lattice. diff --git a/doc/src/fix_rigid_meso.txt b/doc/src/fix_rigid_meso.txt new file mode 100644 index 0000000000..79b014e671 --- /dev/null +++ b/doc/src/fix_rigid_meso.txt @@ -0,0 +1,349 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +fix rigid/meso command :h3 + +[Syntax:] + +fix ID group-ID rigid/meso bodystyle args keyword values ... :pre + +ID, group-ID are documented in "fix"_fix.html command :ulb,l +rigid/meso = style name of this fix command :l +bodystyle = {single} or {molecule} or {group} :l + {single} args = none + {molecule} args = none + {custom} args = {i_propname} or {v_varname} + i_propname = an integer property defined via fix property/atom + v_varname = an atom-style or atomfile-style variable + {group} args = N groupID1 groupID2 ... + N = # of groups + groupID1, groupID2, ... = list of N group IDs :pre + +zero or more keyword/value pairs may be appended :l +keyword = {reinit} or {force} or {torque} or {infile} :l + {reinit} = {yes} or {no} + {force} values = M xflag yflag zflag + M = which rigid body from 1-Nbody (see asterisk form below) + xflag,yflag,zflag = off/on if component of center-of-mass force is active + {torque} values = M xflag yflag zflag + M = which rigid body from 1-Nbody (see asterisk form below) + xflag,yflag,zflag = off/on if component of center-of-mass torque is active + {infile} filename + filename = file with per-body values of mass, center-of-mass, moments of inertia :pre +:ule + +[Examples:] + +fix 1 ellipsoid rigid/meso single +fix 1 rods rigid/meso molecule +fix 1 spheres rigid/meso single force 1 off off on +fix 1 particles rigid/meso molecule force 1*5 off off off force 6*10 off off on +fix 2 spheres rigid/meso group 3 sphere1 sphere2 sphere3 torque * off off off :pre + +[Description:] + +Treat one or more sets of mesoscopic SPH/SDPD particles as independent +rigid bodies. This means that each timestep the total force and torque +on each rigid body is computed as the sum of the forces and torques on +its constituent particles. The coordinates and velocities of the +particles in each body are then updated so that the body moves and +rotates as a single entity using the methods described in the paper by +"(Miller)"_#Miller. Density and internal energy of the particles will +also be updated. This is implemented by creating internal data structures +for each rigid body and performing time integration on these data +structures. Positions and velocities of the constituent particles are +regenerated from the rigid body data structures in every time step. This +restricts which operations and fixes can be applied to rigid bodies. See +below for a detailed discussion. + +The operation of this fix is exactly like that described by the +"fix rigid/nve"_fix_rigid.html command, except that particles' density, +internal energy and extrapolated velocity are also updated. + +NOTE: You should not update the particles in rigid bodies via other +time-integration fixes (e.g. "fix meso"_fix_meso.html, +"fix meso/stationary"_fix_meso_stationary.html), or you will have conflicting +updates to positions and velocities resulting in unphysical behavior in most +cases. When performing a hybrid simulation with some atoms in rigid bodies, +and some not, a separate time integration fix like "fix meso"_fix_meso.html +should be used for the non-rigid particles. + +NOTE: These fixes are overkill if you simply want to hold a collection +of particles stationary or have them move with a constant velocity. To +hold particles stationary use "fix +meso/stationary"_fix_meso_stationary.html instead. If you would like to +move particles with a constant velocity use "fix +meso/move"_fix_meso_move.html. + +IMPORTANT NOTE: The aggregate properties of each rigid body are +calculated at the start of a simulation run and are maintained in +internal data structures. The properties include the position and +velocity of the center-of-mass of the body, its moments of inertia, and +its angular momentum. This is done using the properties of the +constituent particles of the body at that point in time (or see the {infile} +keyword option). Thereafter, changing these properties of individual +particles in the body will have no effect on a rigid body's dynamics, unless +they effect any computation of per-particle forces or torques. If the +keyword {reinit} is set to {yes} (the default), the rigid body data +structures will be recreated at the beginning of each {run} command; +if the keyword {reinit} is set to {no}, the rigid body data structures +will be built only at the very first {run} command and maintained for +as long as the rigid fix is defined. For example, you might think you +could displace the particles in a body or add a large velocity to each particle +in a body to make it move in a desired direction before a 2nd run is +performed, using the "set"_set.html or +"displace_atoms"_displace_atoms.html or "velocity"_velocity.html +commands. But these commands will not affect the internal attributes +of the body unless {reinit} is set to {yes}. With {reinit} set to {no} +(or using the {infile} option, which implies {reinit} {no}) the position +and velocity of individual particles in the body will be reset when time +integration starts again. + +:line + +Each rigid body must have two or more particles. A particle can belong +to at most one rigid body. Which particles are in which bodies can be +defined via several options. + +For bodystyle {single} the entire fix group of particles is treated as +one rigid body. + +For bodystyle {molecule}, particles are grouped into rigid bodies by their +respective molecule IDs: each set of particles in the fix group with the +same molecule ID is treated as a different rigid body. Note that particles +with a molecule ID = 0 will be treated as a single rigid body. For a +system with solvent (typically this is particles with molecule ID = 0) +surrounding rigid bodies, this may not be what you want. Thus you +should be careful to use a fix group that only includes particles you +want to be part of rigid bodies. + +Bodystyle {custom} is similar to bodystyle {molecule} except that it +is more flexible in using other per-atom properties to define the sets +of particles that form rigid bodies. An integer vector defined by the +"fix property/atom"_fix_property_atom.html command can be used. Or an +"atom-style or atomfile-style variable"_variable.html can be used; the +floating-point value produced by the variable is rounded to an +integer. As with bondstyle {molecule}, each set of particles in the fix +groups with the same integer value is treated as a different rigid +body. Since fix property/atom vectors and atom-style variables +produce values for all particles, you should be careful to use a fix group +that only includes particles you want to be part of rigid bodies. + +For bodystyle {group}, each of the listed groups is treated as a +separate rigid body. Only particles that are also in the fix group are +included in each rigid body. + +NOTE: To compute the initial center-of-mass position and other +properties of each rigid body, the image flags for each particle in the +body are used to "unwrap" the particle coordinates. Thus you must +insure that these image flags are consistent so that the unwrapping +creates a valid rigid body (one where the particles are close together) +, particularly if the particles in a single rigid body straddle a +periodic boundary. This means the input data file or restart file must +define the image flags for each particle consistently or that you have +used the "set"_set.html command to specify them correctly. If a +dimension is non-periodic then the image flag of each particle must be +0 in that dimension, else an error is generated. + +By default, each rigid body is acted on by other particles which induce +an external force and torque on its center of mass, causing it to +translate and rotate. Components of the external center-of-mass force +and torque can be turned off by the {force} and {torque} keywords. +This may be useful if you wish a body to rotate but not translate, or +vice versa, or if you wish it to rotate or translate continuously +unaffected by interactions with other particles. Note that if you +expect a rigid body not to move or rotate by using these keywords, you +must insure its initial center-of-mass translational or angular +velocity is 0.0. Otherwise the initial translational or angular +momentum the body has will persist. + +An xflag, yflag, or zflag set to {off} means turn off the component of +force or torque in that dimension. A setting of {on} means turn on +the component, which is the default. Which rigid body(s) the settings +apply to is determined by the first argument of the {force} and +{torque} keywords. It can be an integer M from 1 to Nbody, where +Nbody is the number of rigid bodies defined. A wild-card asterisk can +be used in place of, or in conjunction with, the M argument to set the +flags for multiple rigid bodies. This takes the form "*" or "*n" or +"n*" or "m*n". If N = the number of rigid bodies, then an asterisk +with no numeric values means all bodies from 1 to N. A leading +asterisk means all bodies from 1 to n (inclusive). A trailing +asterisk means all bodies from n to N (inclusive). A middle asterisk +means all bodies from m to n (inclusive). Note that you can use the +{force} or {torque} keywords as many times as you like. If a +particular rigid body has its component flags set multiple times, the +settings from the final keyword are used. + +For computational efficiency, you should typically define one fix +rigid/meso command which includes all the desired rigid bodies. LAMMPS +will allow multiple rigid/meso fixes to be defined, but it is more +expensive. + +:line + +The keyword/value option pairs are used in the following ways. + +The {reinit} keyword determines, whether the rigid body properties +are re-initialized between run commands. With the option {yes} (the +default) this is done, with the option {no} this is not done. Turning +off the re-initialization can be helpful to protect rigid bodies against +unphysical manipulations between runs or when properties cannot be +easily re-computed (e.g. when read from a file). When using the {infile} +keyword, the {reinit} option is automatically set to {no}. + +:line + +The {infile} keyword allows a file of rigid body attributes to be read +in from a file, rather then having LAMMPS compute them. There are 5 +such attributes: the total mass of the rigid body, its center-of-mass +position, its 6 moments of inertia, its center-of-mass velocity, and +the 3 image flags of the center-of-mass position. For rigid bodies +consisting of point particles or non-overlapping finite-size +particles, LAMMPS can compute these values accurately. However, for +rigid bodies consisting of finite-size particles which overlap each +other, LAMMPS will ignore the overlaps when computing these 4 +attributes. The amount of error this induces depends on the amount of +overlap. To avoid this issue, the values can be pre-computed +(e.g. using Monte Carlo integration). + +The format of the file is as follows. Note that the file does not +have to list attributes for every rigid body integrated by fix rigid. +Only bodies which the file specifies will have their computed +attributes overridden. The file can contain initial blank lines or +comment lines starting with "#" which are ignored. The first +non-blank, non-comment line should list N = the number of lines to +follow. The N successive lines contain the following information: + +ID1 masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz vxcm vycm vzcm lx ly lz ixcm iycm izcm +ID2 masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz vxcm vycm vzcm lx ly lz ixcm iycm izcm +... +IDN masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz vxcm vycm vzcm lx ly lz ixcm iycm izcm :pre + +The rigid body IDs are all positive integers. For the {single} +bodystyle, only an ID of 1 can be used. For the {group} bodystyle, +IDs from 1 to Ng can be used where Ng is the number of specified +groups. For the {molecule} bodystyle, use the molecule ID for the +atoms in a specific rigid body as the rigid body ID. + +The masstotal and center-of-mass coordinates (xcm,ycm,zcm) are +self-explanatory. The center-of-mass should be consistent with what +is calculated for the position of the rigid body with all its atoms +unwrapped by their respective image flags. If this produces a +center-of-mass that is outside the simulation box, LAMMPS wraps it +back into the box. + +The 6 moments of inertia (ixx,iyy,izz,ixy,ixz,iyz) should be the +values consistent with the current orientation of the rigid body +around its center of mass. The values are with respect to the +simulation box XYZ axes, not with respect to the principal axes of the +rigid body itself. LAMMPS performs the latter calculation internally. + +The (vxcm,vycm,vzcm) values are the velocity of the center of mass. +The (lx,ly,lz) values are the angular momentum of the body. The +(vxcm,vycm,vzcm) and (lx,ly,lz) values can simply be set to 0 if you +wish the body to have no initial motion. + +The (ixcm,iycm,izcm) values are the image flags of the center of mass +of the body. For periodic dimensions, they specify which image of the +simulation box the body is considered to be in. An image of 0 means +it is inside the box as defined. A value of 2 means add 2 box lengths +to get the true value. A value of -1 means subtract 1 box length to +get the true value. LAMMPS updates these flags as the rigid bodies +cross periodic boundaries during the simulation. + +NOTE: If you use the {infile} keyword and write restart +files during a simulation, then each time a restart file is written, +the fix also write an auxiliary restart file with the name +rfile.rigid, where "rfile" is the name of the restart file, +e.g. tmp.restart.10000 and tmp.restart.10000.rigid. This auxiliary +file is in the same format described above. Thus it can be used in a +new input script that restarts the run and re-specifies a rigid fix +using an {infile} keyword and the appropriate filename. Note that the +auxiliary file will contain one line for every rigid body, even if the +original file only listed a subset of the rigid bodies. + +:line + +[Restart, fix_modify, output, run start/stop, minimize info:] + +No information is written to "binary restart files"_restart.html. +If the {infile} keyword is used, an auxiliary file is written out +with rigid body information each time a restart file is written, as +explained above for the {infile} keyword. + +None of the "fix_modify"_fix_modify.html options are relevant to this +fix. + +This fix computes a global array of values which can be accessed by +various "output commands"_Howto_output.html. + +The number of rows in the array is equal to the number of rigid +bodies. The number of columns is 28. Thus for each rigid body, 28 +values are stored: the xyz coords of the center of mass (COM), the xyz +components of the COM velocity, the xyz components of the force acting +on the COM, the components of the 4-vector quaternion representing the +orientation of the rigid body, the xyz components of the angular momentum +of the body around its COM, the xyz components of the torque acting on the +COM, the 3 principal components of the moment of inertia and the xyz image +flags of the COM. + +The center of mass (COM) for each body is similar to unwrapped +coordinates written to a dump file. It will always be inside (or +slightly outside) the simulation box. The image flags have the same +meaning as image flags for particle positions (see the "dump" command). +This means you can calculate the unwrapped COM by applying the image +flags to the COM, the same as when unwrapped coordinates are written +to a dump file. + +The force and torque values in the array are not affected by the +{force} and {torque} keywords in the fix rigid command; they reflect +values before any changes are made by those keywords. + +The ordering of the rigid bodies (by row in the array) is as follows. +For the {single} keyword there is just one rigid body. For the +{molecule} keyword, the bodies are ordered by ascending molecule ID. +For the {group} keyword, the list of group IDs determines the ordering +of bodies. + +The array values calculated by this fix are "intensive", meaning they +are independent of the number of particles in the simulation. + +No parameter of this fix can be used with the {start/stop} keywords of +the "run"_run.html command. + +This fix is not invoked during "energy minimization"_minimize.html. + +:line + +[Restrictions:] + +This fix is part of the USER-SDPD package and also depends on the RIGID +package. It is only enabled if LAMMPS was built with both packages. See +the "Build package"_Build_package.html doc page for more info. + +This fix requires that atoms store density and internal energy as +defined by the "atom_style meso"_atom_style.html command. + +All particles in the group must be mesoscopic SPH/SDPD particles. + +[Related commands:] + +"fix meso/move"_fix_meso_move.html, "fix rigid"_fix_rigid.html, +"neigh_modify exclude"_neigh_modify.html + +[Default:] + +The option defaults are force * on on on and torque * on on on, +meaning all rigid bodies are acted on by center-of-mass force and +torque. Also reinit = yes. + +:line + +:link(Miller) +[(Miller)] Miller, Eleftheriou, Pattnaik, Ndirango, and Newns, +J Chem Phys, 116, 8649 (2002). diff --git a/doc/src/fixes.txt b/doc/src/fixes.txt index 66f6633124..36ccf4269a 100644 --- a/doc/src/fixes.txt +++ b/doc/src/fixes.txt @@ -73,6 +73,7 @@ Fixes :h1 fix_lineforce fix_manifoldforce fix_meso + fix_meso_move fix_meso_stationary fix_momentum fix_move @@ -137,6 +138,7 @@ Fixes :h1 fix_restrain fix_rhok fix_rigid + fix_rigid_meso fix_rx fix_saed_vtk fix_setforce diff --git a/doc/src/lammps.book b/doc/src/lammps.book index f739712f08..cb328c4089 100644 --- a/doc/src/lammps.book +++ b/doc/src/lammps.book @@ -293,6 +293,7 @@ fix_lb_viscous.html fix_lineforce.html fix_manifoldforce.html fix_meso.html +fix_meso_move.html fix_meso_stationary.html fix_momentum.html fix_move.html @@ -356,6 +357,7 @@ fix_reaxc_species.html fix_recenter.html fix_restrain.html fix_rigid.html +fix_rigid_meso.html fix_rhok.html fix_rx.html fix_saed_vtk.html @@ -615,6 +617,7 @@ pair_reax.html pair_reaxc.html pair_resquared.html pair_sdk.html +pair_sdpd_taitwater_isothermal.html pair_smd_hertz.html pair_smd_tlsph.html pair_smd_triangulated_surface.html diff --git a/doc/src/pair_sdpd_taitwater_isothermal.txt b/doc/src/pair_sdpd_taitwater_isothermal.txt new file mode 100644 index 0000000000..97b68ceba7 --- /dev/null +++ b/doc/src/pair_sdpd_taitwater_isothermal.txt @@ -0,0 +1,108 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +pair_style sdpd/taitwater/isothermal command :h3 + +[Syntax:] + +pair_style sdpd/taitwater/isothermal temperature viscosity seed +:pre + +temperature = temperature of the fluid (temperature units) +viscosity = dynamic viscosity of the fluid (mass*distance/time units) +seed = random number generator seed (positive integer, optional) :ul + +[Examples:] + +pair_style sdpd/taitwater/isothermal 300. 1. 28681 +pair_coeff * * 1000.0 1430.0 2.4 :pre + +[Description:] + +The sdpd/taitwater/isothermal style computes forces between mesoscopic +particles according to the Smoothed Dissipative Particle Dynamics model +described in this paper by "(Español and Revenga)"_#Español_Revenga under +the following assumptions: + +:olb +The temperature is constant and uniform. :l +The shear viscosity is constant and uniform. :l +The volume viscosity is negligible before the shear viscosity. :l +The Boltzmann constant is negligible before the heat capacity of a +single mesoscopic particle of fluid. :ole,l + +The third assumption is true for water in nearly incompressible flows. +The fourth holds true for water for any reasonable size one can +imagine for a mesoscopic particle. + +The pressure forces between particles will be computed according to +Tait's equation of state: + +:c,image(Eqs/pair_sph_tait.jpg) + +where gamma = 7 and B = c_0^2 rho_0 / gamma, with rho_0 being the +reference density and c_0 the reference speed of sound. + +The laminar viscosity and the random forces will be computed according +to formulas described in "(Español and Revenga)"_#Español_Revenga. + +IMPORTANT NOTE: Similar to "brownian"_pair_brownian.html and +"dpd"_pair_dpd.html styles, the "newton"_newton.html setting for +pairwise interactions needs to be on when running LAMMPS in parallel +if you want to ensure linear momentum conservation. Otherwise random +forces generated for pairs straddling processor boundary will not be +equal and opposite. + +NOTE: The actual random seed used will be a mix of what you specify +and other parameters like the MPI ranks. This is to ensure that +different MPI tasks have distinct seeds. + +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. + +rho0 reference density (mass/volume units) +c0 reference soundspeed (distance/time units) +h kernel function cutoff (distance units) :ul + +:line + +[Mixing, shift, table, tail correction, restart, rRESPA info]: + +This style does not support mixing. Thus, coefficients for all +I,J pairs must be specified explicitly. + +This style does not support the "pair_modify"_pair_modify.html +shift, table, and tail options. + +This style does not write information to "binary restart +files"_restart.html. Thus, you need to re-specify the pair_style and +pair_coeff commands in an input script that reads a restart file. + +This style can only be used via the {pair} keyword of the "run_style +respa"_run_style.html command. It does not support the {inner}, +{middle}, {outer} keywords. + +[Restrictions:] + +This pair style is part of the USER-SDPD package. It is only enabled +if LAMMPS was built with that package. See the "Build +package"_Build_package.html doc page for more info. + +[Related commands:] + +"pair coeff"_pair_coeff.html, "pair sph/rhosum"_pair_sph_rhosum.html + +[Default:] + +The default seed is 0 (before mixing). + +:line + +:link(Español_Revenga) +[(Español and Revenga)] Español, Revenga, Physical Review E, 67, 026705 (2003). diff --git a/doc/src/pair_style.txt b/doc/src/pair_style.txt index 14e3f79215..f6098fa005 100644 --- a/doc/src/pair_style.txt +++ b/doc/src/pair_style.txt @@ -268,6 +268,7 @@ pair"_Commands_pair.html doc page are followed by one or more of "reax/c"_pair_reaxc.html - ReaxFF potential in C "rebo"_pair_airebo.html - 2nd generation REBO potential of Brenner "resquared"_pair_resquared.html - Everaers RE-Squared ellipsoidal potential +"sdpd/taitwater/isothermal"_pair_sdpd_taitwater_isothermal.html - smoothed dissipative particle dynamics for water at isothermal conditions "smd/hertz"_pair_smd_hertz.html - "smd/tlsph"_pair_smd_tlsph.html - "smd/tri_surface"_pair_smd_triangulated_surface.html - diff --git a/doc/src/pairs.txt b/doc/src/pairs.txt index d535798482..ca79051053 100644 --- a/doc/src/pairs.txt +++ b/doc/src/pairs.txt @@ -86,6 +86,7 @@ Pair Styles :h1 pair_reaxc pair_resquared pair_sdk + pair_sdpd_taitwater_isothermal pair_smd_hertz pair_smd_tlsph pair_smd_triangulated_surface diff --git a/examples/USER/sdpd/2d-diffusion-in-shear-flow/in.lammps b/examples/USER/sdpd/2d-diffusion-in-shear-flow/in.lammps new file mode 100644 index 0000000000..138f7b5338 --- /dev/null +++ b/examples/USER/sdpd/2d-diffusion-in-shear-flow/in.lammps @@ -0,0 +1,61 @@ +dimension 2 +units micro +atom_style meso + +variable R equal 0.5 # radius of sphere micrometers +variable a equal $R/5 # lattice spacing micrometers +variable Lf equal $R*3 +variable Lb equal $R*4 +variable wall_velocity equal 0.01 # micrometers/microsecond +variable T equal 300. +variable rho_0 equal 1. # density picograms/micrometer^3 +variable c_0 equal 100. # speed of sound micrometers/microsecond +variable mu equal 1. # dynamic viscosity picogram/(micrometer-microsecond) +variable h equal $a*4.5 # kernel function cutoff micrometers +variable mass equal $a*$a*$a*${rho_0} +variable dt equal 1e-3 # timestep microseconds +variable skin equal 0.2*$h + +region box block -${Lb} ${Lb} -${Lb} ${Lb} 0 ${a} units box +create_box 4 box +lattice sq $a + +create_atoms 1 box + +region sphere sphere 0 0 0 $R units box +set region sphere type 2 + +region upper_wall block INF INF +${Lf} INF INF INF units box +set region upper_wall type 3 + +region lower_wall block INF INF INF -${Lf} INF INF units box +set region lower_wall type 4 + +group fluid type 1 +group sphere type 2 +group upper_wall type 3 +group lower_wall type 4 + +mass * ${mass} +set group all meso/rho ${rho_0} + +pair_style sdpd/taitwater/isothermal $T ${mu} 76787 # temperature viscosity random_seed +pair_coeff * * ${rho_0} ${c_0} ${h} + +fix 1 fluid meso +fix 2 sphere rigid/meso single +fix 3 upper_wall meso/move linear +${wall_velocity} 0 0 units box +fix 4 lower_wall meso/move linear -${wall_velocity} 0 0 units box + +fix 2d all enforce2d + +neighbor ${skin} bin +neigh_modify delay 0 every 1 check yes +timestep ${dt} + +dump dump_id all atom 100 dump.lammpstrj + +thermo 100 +thermo_style custom step time nbuild ndanger + +run 10000 diff --git a/examples/USER/sdpd/2d-diffusion-in-shear-flow/log.24Oct18.2d-diffusion-in-shear-flow.g++.1 b/examples/USER/sdpd/2d-diffusion-in-shear-flow/log.24Oct18.2d-diffusion-in-shear-flow.g++.1 new file mode 100644 index 0000000000..cd01601292 --- /dev/null +++ b/examples/USER/sdpd/2d-diffusion-in-shear-flow/log.24Oct18.2d-diffusion-in-shear-flow.g++.1 @@ -0,0 +1,247 @@ +LAMMPS (24 Oct 2018) +dimension 2 +units micro +atom_style meso + +variable R equal 0.5 # radius of sphere micrometers +variable a equal $R/5 # lattice spacing micrometers +variable a equal 0.5/5 +variable Lf equal $R*3 +variable Lf equal 0.5*3 +variable Lb equal $R*4 +variable Lb equal 0.5*4 +variable wall_velocity equal 0.01 # micrometers/microsecond +variable T equal 300. +variable rho_0 equal 1. # density picograms/micrometer^3 +variable c_0 equal 100. # speed of sound micrometers/microsecond +variable mu equal 1. # dynamic viscosity picogram/(micrometer-microsecond) +variable h equal $a*4.5 # kernel function cutoff micrometers +variable h equal 0.1*4.5 +variable mass equal $a*$a*$a*${rho_0} +variable mass equal 0.1*$a*$a*${rho_0} +variable mass equal 0.1*0.1*$a*${rho_0} +variable mass equal 0.1*0.1*0.1*${rho_0} +variable mass equal 0.1*0.1*0.1*1 +variable dt equal 1e-3 # timestep microseconds +variable skin equal 0.2*$h +variable skin equal 0.2*0.45 + +region box block -${Lb} ${Lb} -${Lb} ${Lb} 0 ${a} units box +region box block -2 ${Lb} -${Lb} ${Lb} 0 ${a} units box +region box block -2 2 -${Lb} ${Lb} 0 ${a} units box +region box block -2 2 -2 ${Lb} 0 ${a} units box +region box block -2 2 -2 2 0 ${a} units box +region box block -2 2 -2 2 0 0.1 units box +create_box 4 box +Created orthogonal box = (-2 -2 0) to (2 2 0.1) + 1 by 1 by 1 MPI processor grid +lattice sq $a +lattice sq 0.1 +Lattice spacing in x,y,z = 0.1 0.1 0.1 + +create_atoms 1 box +Created 1600 atoms + Time spent = 0.00169706 secs + +region sphere sphere 0 0 0 $R units box +region sphere sphere 0 0 0 0.5 units box +set region sphere type 2 + 81 settings made for type + +region upper_wall block INF INF +${Lf} INF INF INF units box +region upper_wall block INF INF +1.5 INF INF INF units box +set region upper_wall type 3 + 200 settings made for type + +region lower_wall block INF INF INF -${Lf} INF INF units box +region lower_wall block INF INF INF -1.5 INF INF units box +set region lower_wall type 4 + 240 settings made for type + +group fluid type 1 +1079 atoms in group fluid +group sphere type 2 +81 atoms in group sphere +group upper_wall type 3 +200 atoms in group upper_wall +group lower_wall type 4 +240 atoms in group lower_wall + +mass * ${mass} +mass * 0.001 +set group all meso/rho ${rho_0} +set group all meso/rho 1 + 1600 settings made for meso/rho + +pair_style sdpd/taitwater/isothermal $T ${mu} 76787 # temperature viscosity random_seed +pair_style sdpd/taitwater/isothermal 300 ${mu} 76787 +pair_style sdpd/taitwater/isothermal 300 1 76787 +pair_coeff * * ${rho_0} ${c_0} ${h} +pair_coeff * * 1 ${c_0} ${h} +pair_coeff * * 1 100 ${h} +pair_coeff * * 1 100 0.45 + +fix 1 fluid meso +fix 2 sphere rigid/meso single +1 rigid bodies with 81 atoms +fix 3 upper_wall meso/move linear +${wall_velocity} 0 0 units box +fix 3 upper_wall meso/move linear +0.01 0 0 units box +fix 4 lower_wall meso/move linear -${wall_velocity} 0 0 units box +fix 4 lower_wall meso/move linear -0.01 0 0 units box + +fix 2d all enforce2d + +neighbor ${skin} bin +neighbor 0.09 bin +neigh_modify delay 0 every 1 check yes +timestep ${dt} +timestep 0.001 + +dump dump_id all atom 100 dump.lammpstrj + +thermo 100 +thermo_style custom step time nbuild ndanger + +run 10000 +Neighbor list info ... + update every 1 steps, delay 0 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 0.54 + ghost atom cutoff = 0.54 + binsize = 0.27, bins = 15 15 1 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair sdpd/taitwater/isothermal, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/2d/newton + bin: standard +Per MPI rank memory allocation (min/avg/max) = 6.937 | 6.937 | 6.937 Mbytes +Step Time Nbuild Ndanger + 0 0 0 0 + 100 0.1 0 0 + 200 0.2 0 0 + 300 0.3 0 0 + 400 0.4 0 0 + 500 0.5 1 0 + 600 0.6 1 0 + 700 0.7 2 0 + 800 0.8 2 0 + 900 0.9 2 0 + 1000 1 3 0 + 1100 1.1 3 0 + 1200 1.2 3 0 + 1300 1.3 4 0 + 1400 1.4 4 0 + 1500 1.5 4 0 + 1600 1.6 5 0 + 1700 1.7 5 0 + 1800 1.8 5 0 + 1900 1.9 6 0 + 2000 2 6 0 + 2100 2.1 6 0 + 2200 2.2 7 0 + 2300 2.3 7 0 + 2400 2.4 7 0 + 2500 2.5 8 0 + 2600 2.6 8 0 + 2700 2.7 8 0 + 2800 2.8 9 0 + 2900 2.9 9 0 + 3000 3 9 0 + 3100 3.1 10 0 + 3200 3.2 10 0 + 3300 3.3 10 0 + 3400 3.4 11 0 + 3500 3.5 11 0 + 3600 3.6 11 0 + 3700 3.7 12 0 + 3800 3.8 12 0 + 3900 3.9 12 0 + 4000 4 13 0 + 4100 4.1 13 0 + 4200 4.2 14 0 + 4300 4.3 14 0 + 4400 4.4 14 0 + 4500 4.5 15 0 + 4600 4.6 15 0 + 4700 4.7 15 0 + 4800 4.8 16 0 + 4900 4.9 16 0 + 5000 5 16 0 + 5100 5.1 17 0 + 5200 5.2 17 0 + 5300 5.3 17 0 + 5400 5.4 17 0 + 5500 5.5 18 0 + 5600 5.6 18 0 + 5700 5.7 19 0 + 5800 5.8 19 0 + 5900 5.9 19 0 + 6000 6 20 0 + 6100 6.1 20 0 + 6200 6.2 21 0 + 6300 6.3 21 0 + 6400 6.4 21 0 + 6500 6.5 22 0 + 6600 6.6 22 0 + 6700 6.7 22 0 + 6800 6.8 23 0 + 6900 6.9 23 0 + 7000 7 23 0 + 7100 7.1 24 0 + 7200 7.2 24 0 + 7300 7.3 25 0 + 7400 7.4 25 0 + 7500 7.5 25 0 + 7600 7.6 26 0 + 7700 7.7 26 0 + 7800 7.8 26 0 + 7900 7.9 27 0 + 8000 8 27 0 + 8100 8.1 27 0 + 8200 8.2 28 0 + 8300 8.3 28 0 + 8400 8.4 28 0 + 8500 8.5 29 0 + 8600 8.6 29 0 + 8700 8.7 30 0 + 8800 8.8 30 0 + 8900 8.9 30 0 + 9000 9 31 0 + 9100 9.1 31 0 + 9200 9.2 31 0 + 9300 9.3 32 0 + 9400 9.4 32 0 + 9500 9.5 32 0 + 9600 9.6 33 0 + 9700 9.7 33 0 + 9800 9.8 33 0 + 9900 9.9 34 0 + 10000 10 34 0 +Loop time of 144.208 on 1 procs for 10000 steps with 1600 atoms + +Performance: 5991348.580 ns/day, 0.000 hours/ns, 69.344 timesteps/s +99.7% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 143.08 | 143.08 | 143.08 | 0.0 | 99.22 +Neigh | 0.033195 | 0.033195 | 0.033195 | 0.0 | 0.02 +Comm | 0.24139 | 0.24139 | 0.24139 | 0.0 | 0.17 +Output | 0.11687 | 0.11687 | 0.11687 | 0.0 | 0.08 +Modify | 0.61566 | 0.61566 | 0.61566 | 0.0 | 0.43 +Other | | 0.117 | | | 0.08 + +Nlocal: 1600 ave 1600 max 1600 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 993 ave 993 max 993 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 73236 ave 73236 max 73236 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 73236 +Ave neighs/atom = 45.7725 +Neighbor list builds = 34 +Dangerous builds = 0 +Total wall time: 0:02:24 diff --git a/examples/USER/sdpd/2d-diffusion-in-shear-flow/log.24Oct18.2d-diffusion-in-shear-flow.g++.4 b/examples/USER/sdpd/2d-diffusion-in-shear-flow/log.24Oct18.2d-diffusion-in-shear-flow.g++.4 new file mode 100644 index 0000000000..77823e00fc --- /dev/null +++ b/examples/USER/sdpd/2d-diffusion-in-shear-flow/log.24Oct18.2d-diffusion-in-shear-flow.g++.4 @@ -0,0 +1,247 @@ +LAMMPS (24 Oct 2018) +dimension 2 +units micro +atom_style meso + +variable R equal 0.5 # radius of sphere micrometers +variable a equal $R/5 # lattice spacing micrometers +variable a equal 0.5/5 +variable Lf equal $R*3 +variable Lf equal 0.5*3 +variable Lb equal $R*4 +variable Lb equal 0.5*4 +variable wall_velocity equal 0.01 # micrometers/microsecond +variable T equal 300. +variable rho_0 equal 1. # density picograms/micrometer^3 +variable c_0 equal 100. # speed of sound micrometers/microsecond +variable mu equal 1. # dynamic viscosity picogram/(micrometer-microsecond) +variable h equal $a*4.5 # kernel function cutoff micrometers +variable h equal 0.1*4.5 +variable mass equal $a*$a*$a*${rho_0} +variable mass equal 0.1*$a*$a*${rho_0} +variable mass equal 0.1*0.1*$a*${rho_0} +variable mass equal 0.1*0.1*0.1*${rho_0} +variable mass equal 0.1*0.1*0.1*1 +variable dt equal 1e-3 # timestep microseconds +variable skin equal 0.2*$h +variable skin equal 0.2*0.45 + +region box block -${Lb} ${Lb} -${Lb} ${Lb} 0 ${a} units box +region box block -2 ${Lb} -${Lb} ${Lb} 0 ${a} units box +region box block -2 2 -${Lb} ${Lb} 0 ${a} units box +region box block -2 2 -2 ${Lb} 0 ${a} units box +region box block -2 2 -2 2 0 ${a} units box +region box block -2 2 -2 2 0 0.1 units box +create_box 4 box +Created orthogonal box = (-2 -2 0) to (2 2 0.1) + 2 by 2 by 1 MPI processor grid +lattice sq $a +lattice sq 0.1 +Lattice spacing in x,y,z = 0.1 0.1 0.1 + +create_atoms 1 box +Created 1600 atoms + Time spent = 0.000589566 secs + +region sphere sphere 0 0 0 $R units box +region sphere sphere 0 0 0 0.5 units box +set region sphere type 2 + 81 settings made for type + +region upper_wall block INF INF +${Lf} INF INF INF units box +region upper_wall block INF INF +1.5 INF INF INF units box +set region upper_wall type 3 + 200 settings made for type + +region lower_wall block INF INF INF -${Lf} INF INF units box +region lower_wall block INF INF INF -1.5 INF INF units box +set region lower_wall type 4 + 240 settings made for type + +group fluid type 1 +1079 atoms in group fluid +group sphere type 2 +81 atoms in group sphere +group upper_wall type 3 +200 atoms in group upper_wall +group lower_wall type 4 +240 atoms in group lower_wall + +mass * ${mass} +mass * 0.001 +set group all meso/rho ${rho_0} +set group all meso/rho 1 + 1600 settings made for meso/rho + +pair_style sdpd/taitwater/isothermal $T ${mu} 76787 # temperature viscosity random_seed +pair_style sdpd/taitwater/isothermal 300 ${mu} 76787 +pair_style sdpd/taitwater/isothermal 300 1 76787 +pair_coeff * * ${rho_0} ${c_0} ${h} +pair_coeff * * 1 ${c_0} ${h} +pair_coeff * * 1 100 ${h} +pair_coeff * * 1 100 0.45 + +fix 1 fluid meso +fix 2 sphere rigid/meso single +1 rigid bodies with 81 atoms +fix 3 upper_wall meso/move linear +${wall_velocity} 0 0 units box +fix 3 upper_wall meso/move linear +0.01 0 0 units box +fix 4 lower_wall meso/move linear -${wall_velocity} 0 0 units box +fix 4 lower_wall meso/move linear -0.01 0 0 units box + +fix 2d all enforce2d + +neighbor ${skin} bin +neighbor 0.09 bin +neigh_modify delay 0 every 1 check yes +timestep ${dt} +timestep 0.001 + +dump dump_id all atom 100 dump.lammpstrj + +thermo 100 +thermo_style custom step time nbuild ndanger + +run 10000 +Neighbor list info ... + update every 1 steps, delay 0 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 0.54 + ghost atom cutoff = 0.54 + binsize = 0.27, bins = 15 15 1 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair sdpd/taitwater/isothermal, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/2d/newton + bin: standard +Per MPI rank memory allocation (min/avg/max) = 6.854 | 6.854 | 6.854 Mbytes +Step Time Nbuild Ndanger + 0 0 0 0 + 100 0.1 0 0 + 200 0.2 0 0 + 300 0.3 0 0 + 400 0.4 1 0 + 500 0.5 1 0 + 600 0.6 1 0 + 700 0.7 2 0 + 800 0.8 2 0 + 900 0.9 2 0 + 1000 1 3 0 + 1100 1.1 3 0 + 1200 1.2 4 0 + 1300 1.3 4 0 + 1400 1.4 4 0 + 1500 1.5 4 0 + 1600 1.6 5 0 + 1700 1.7 5 0 + 1800 1.8 5 0 + 1900 1.9 6 0 + 2000 2 6 0 + 2100 2.1 6 0 + 2200 2.2 6 0 + 2300 2.3 7 0 + 2400 2.4 7 0 + 2500 2.5 7 0 + 2600 2.6 8 0 + 2700 2.7 8 0 + 2800 2.8 8 0 + 2900 2.9 9 0 + 3000 3 9 0 + 3100 3.1 9 0 + 3200 3.2 10 0 + 3300 3.3 10 0 + 3400 3.4 10 0 + 3500 3.5 11 0 + 3600 3.6 11 0 + 3700 3.7 11 0 + 3800 3.8 12 0 + 3900 3.9 12 0 + 4000 4 12 0 + 4100 4.1 13 0 + 4200 4.2 13 0 + 4300 4.3 13 0 + 4400 4.4 14 0 + 4500 4.5 14 0 + 4600 4.6 15 0 + 4700 4.7 15 0 + 4800 4.8 15 0 + 4900 4.9 16 0 + 5000 5 16 0 + 5100 5.1 17 0 + 5200 5.2 17 0 + 5300 5.3 17 0 + 5400 5.4 17 0 + 5500 5.5 18 0 + 5600 5.6 18 0 + 5700 5.7 18 0 + 5800 5.8 19 0 + 5900 5.9 19 0 + 6000 6 20 0 + 6100 6.1 20 0 + 6200 6.2 20 0 + 6300 6.3 21 0 + 6400 6.4 21 0 + 6500 6.5 21 0 + 6600 6.6 22 0 + 6700 6.7 22 0 + 6800 6.8 22 0 + 6900 6.9 23 0 + 7000 7 23 0 + 7100 7.1 23 0 + 7200 7.2 24 0 + 7300 7.3 24 0 + 7400 7.4 25 0 + 7500 7.5 25 0 + 7600 7.6 25 0 + 7700 7.7 25 0 + 7800 7.8 26 0 + 7900 7.9 26 0 + 8000 8 26 0 + 8100 8.1 27 0 + 8200 8.2 27 0 + 8300 8.3 27 0 + 8400 8.4 28 0 + 8500 8.5 28 0 + 8600 8.6 28 0 + 8700 8.7 29 0 + 8800 8.8 29 0 + 8900 8.9 29 0 + 9000 9 30 0 + 9100 9.1 30 0 + 9200 9.2 31 0 + 9300 9.3 31 0 + 9400 9.4 31 0 + 9500 9.5 32 0 + 9600 9.6 32 0 + 9700 9.7 32 0 + 9800 9.8 32 0 + 9900 9.9 33 0 + 10000 10 33 0 +Loop time of 63.2372 on 4 procs for 10000 steps with 1600 atoms + +Performance: 13662841.706 ns/day, 0.000 hours/ns, 158.135 timesteps/s +94.3% CPU use with 4 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 51.576 | 53.662 | 55.484 | 23.9 | 84.86 +Neigh | 0.011519 | 0.012395 | 0.013405 | 0.7 | 0.02 +Comm | 6.8389 | 8.5423 | 10.517 | 56.1 | 13.51 +Output | 0.12342 | 0.12513 | 0.1302 | 0.8 | 0.20 +Modify | 0.58708 | 0.69128 | 0.78806 | 11.3 | 1.09 +Other | | 0.2038 | | | 0.32 + +Nlocal: 400 ave 411 max 388 min +Histogram: 1 1 0 0 0 0 0 0 0 2 +Nghost: 552.25 ave 567 max 539 min +Histogram: 2 0 0 0 0 0 0 0 1 1 +Neighs: 18298.8 ave 18781 max 17829 min +Histogram: 2 0 0 0 0 0 0 0 0 2 + +Total # of neighbors = 73195 +Ave neighs/atom = 45.7469 +Neighbor list builds = 33 +Dangerous builds = 0 +Total wall time: 0:01:03 diff --git a/examples/USER/sdpd/2d-diffusion/in.lammps b/examples/USER/sdpd/2d-diffusion/in.lammps new file mode 100644 index 0000000000..6ef36a0cf6 --- /dev/null +++ b/examples/USER/sdpd/2d-diffusion/in.lammps @@ -0,0 +1,49 @@ +dimension 2 +units micro +atom_style meso + +variable R equal 0.5 # radius of sphere micrometers +variable a equal $R/5 # lattice spacing micrometers +variable L equal $R*3 +variable T equal 300. +variable rho_0 equal 1. # density picograms/micrometer^3 +variable c_0 equal 100. # speed of sound micrometers/microsecond +variable mu equal 1. # dynamic viscosity picogram/(micrometer-microsecond) +variable h equal $a*4.5 # kernel function cutoff micrometers +variable mass equal $a*$a*$a*${rho_0} +variable dt equal 1e-3 # timestep microseconds +variable skin equal 0.2*$h + +region box block -$L $L -$L $L 0 $a units box +create_box 2 box +lattice sq $a + +create_atoms 1 box + +region sphere sphere 0 0 0 $R units box +set region sphere type 2 + +group fluid type 1 +group sphere type 2 + +mass * ${mass} +set group all meso/rho ${rho_0} + +pair_style sdpd/taitwater/isothermal $T ${mu} 76787 # temperature viscosity random_seed +pair_coeff * * ${rho_0} ${c_0} ${h} + +fix 1 fluid meso +fix 2 sphere rigid/meso single + +fix 2d all enforce2d + +neighbor ${skin} bin +neigh_modify delay 0 every 1 check yes +timestep ${dt} + +dump dump_id all atom 100 dump.lammpstrj + +thermo 100 +thermo_style custom step time nbuild ndanger + +run 10000 diff --git a/examples/USER/sdpd/2d-diffusion/log.24Oct18.2d-diffusion.g++.1 b/examples/USER/sdpd/2d-diffusion/log.24Oct18.2d-diffusion.g++.1 new file mode 100644 index 0000000000..d44c0fd6f4 --- /dev/null +++ b/examples/USER/sdpd/2d-diffusion/log.24Oct18.2d-diffusion.g++.1 @@ -0,0 +1,226 @@ +LAMMPS (24 Oct 2018) +dimension 2 +units micro +atom_style meso + +variable R equal 0.5 # radius of sphere micrometers +variable a equal $R/5 # lattice spacing micrometers +variable a equal 0.5/5 +variable L equal $R*3 +variable L equal 0.5*3 +variable T equal 300. +variable rho_0 equal 1. # density picograms/micrometer^3 +variable c_0 equal 100. # speed of sound micrometers/microsecond +variable mu equal 1. # dynamic viscosity picogram/(micrometer-microsecond) +variable h equal $a*4.5 # kernel function cutoff micrometers +variable h equal 0.1*4.5 +variable mass equal $a*$a*$a*${rho_0} +variable mass equal 0.1*$a*$a*${rho_0} +variable mass equal 0.1*0.1*$a*${rho_0} +variable mass equal 0.1*0.1*0.1*${rho_0} +variable mass equal 0.1*0.1*0.1*1 +variable dt equal 1e-3 # timestep microseconds +variable skin equal 0.2*$h +variable skin equal 0.2*0.45 + +region box block -$L $L -$L $L 0 $a units box +region box block -1.5 $L -$L $L 0 $a units box +region box block -1.5 1.5 -$L $L 0 $a units box +region box block -1.5 1.5 -1.5 $L 0 $a units box +region box block -1.5 1.5 -1.5 1.5 0 $a units box +region box block -1.5 1.5 -1.5 1.5 0 0.1 units box +create_box 2 box +Created orthogonal box = (-1.5 -1.5 0) to (1.5 1.5 0.1) + 1 by 1 by 1 MPI processor grid +lattice sq $a +lattice sq 0.1 +Lattice spacing in x,y,z = 0.1 0.1 0.1 + +create_atoms 1 box +Created 900 atoms + Time spent = 0.0015769 secs + +region sphere sphere 0 0 0 $R units box +region sphere sphere 0 0 0 0.5 units box +set region sphere type 2 + 81 settings made for type + +group fluid type 1 +819 atoms in group fluid +group sphere type 2 +81 atoms in group sphere + +mass * ${mass} +mass * 0.001 +set group all meso/rho ${rho_0} +set group all meso/rho 1 + 900 settings made for meso/rho + +pair_style sdpd/taitwater/isothermal $T ${mu} 76787 # temperature viscosity random_seed +pair_style sdpd/taitwater/isothermal 300 ${mu} 76787 +pair_style sdpd/taitwater/isothermal 300 1 76787 +pair_coeff * * ${rho_0} ${c_0} ${h} +pair_coeff * * 1 ${c_0} ${h} +pair_coeff * * 1 100 ${h} +pair_coeff * * 1 100 0.45 + +fix 1 fluid meso +fix 2 sphere rigid/meso single +1 rigid bodies with 81 atoms + +fix 2d all enforce2d + +neighbor ${skin} bin +neighbor 0.09 bin +neigh_modify delay 0 every 1 check yes +timestep ${dt} +timestep 0.001 + +dump dump_id all atom 100 dump.lammpstrj + +thermo 100 +thermo_style custom step time nbuild ndanger + +run 10000 +Neighbor list info ... + update every 1 steps, delay 0 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 0.54 + ghost atom cutoff = 0.54 + binsize = 0.27, bins = 12 12 1 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair sdpd/taitwater/isothermal, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/2d/newton + bin: standard +Per MPI rank memory allocation (min/avg/max) = 6.137 | 6.137 | 6.137 Mbytes +Step Time Nbuild Ndanger + 0 0 0 0 + 100 0.1 0 0 + 200 0.2 0 0 + 300 0.3 0 0 + 400 0.4 1 0 + 500 0.5 1 0 + 600 0.6 1 0 + 700 0.7 2 0 + 800 0.8 2 0 + 900 0.9 2 0 + 1000 1 3 0 + 1100 1.1 3 0 + 1200 1.2 3 0 + 1300 1.3 4 0 + 1400 1.4 4 0 + 1500 1.5 4 0 + 1600 1.6 5 0 + 1700 1.7 5 0 + 1800 1.8 6 0 + 1900 1.9 6 0 + 2000 2 6 0 + 2100 2.1 7 0 + 2200 2.2 7 0 + 2300 2.3 7 0 + 2400 2.4 7 0 + 2500 2.5 8 0 + 2600 2.6 8 0 + 2700 2.7 8 0 + 2800 2.8 9 0 + 2900 2.9 9 0 + 3000 3 10 0 + 3100 3.1 10 0 + 3200 3.2 10 0 + 3300 3.3 11 0 + 3400 3.4 11 0 + 3500 3.5 11 0 + 3600 3.6 12 0 + 3700 3.7 12 0 + 3800 3.8 12 0 + 3900 3.9 13 0 + 4000 4 13 0 + 4100 4.1 13 0 + 4200 4.2 14 0 + 4300 4.3 14 0 + 4400 4.4 14 0 + 4500 4.5 15 0 + 4600 4.6 15 0 + 4700 4.7 15 0 + 4800 4.8 16 0 + 4900 4.9 16 0 + 5000 5 17 0 + 5100 5.1 17 0 + 5200 5.2 17 0 + 5300 5.3 17 0 + 5400 5.4 18 0 + 5500 5.5 18 0 + 5600 5.6 18 0 + 5700 5.7 19 0 + 5800 5.8 19 0 + 5900 5.9 19 0 + 6000 6 19 0 + 6100 6.1 20 0 + 6200 6.2 20 0 + 6300 6.3 20 0 + 6400 6.4 21 0 + 6500 6.5 21 0 + 6600 6.6 21 0 + 6700 6.7 21 0 + 6800 6.8 22 0 + 6900 6.9 22 0 + 7000 7 22 0 + 7100 7.1 23 0 + 7200 7.2 23 0 + 7300 7.3 23 0 + 7400 7.4 24 0 + 7500 7.5 24 0 + 7600 7.6 24 0 + 7700 7.7 25 0 + 7800 7.8 25 0 + 7900 7.9 26 0 + 8000 8 26 0 + 8100 8.1 26 0 + 8200 8.2 26 0 + 8300 8.3 27 0 + 8400 8.4 27 0 + 8500 8.5 27 0 + 8600 8.6 28 0 + 8700 8.7 28 0 + 8800 8.8 28 0 + 8900 8.9 29 0 + 9000 9 29 0 + 9100 9.1 29 0 + 9200 9.2 30 0 + 9300 9.3 30 0 + 9400 9.4 30 0 + 9500 9.5 30 0 + 9600 9.6 31 0 + 9700 9.7 31 0 + 9800 9.8 32 0 + 9900 9.9 32 0 + 10000 10 32 0 +Loop time of 80.9456 on 1 procs for 10000 steps with 900 atoms + +Performance: 10673829.855 ns/day, 0.000 hours/ns, 123.540 timesteps/s +99.8% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 80.306 | 80.306 | 80.306 | 0.0 | 99.21 +Neigh | 0.017418 | 0.017418 | 0.017418 | 0.0 | 0.02 +Comm | 0.16939 | 0.16939 | 0.16939 | 0.0 | 0.21 +Output | 0.070281 | 0.070281 | 0.070281 | 0.0 | 0.09 +Modify | 0.3154 | 0.3154 | 0.3154 | 0.0 | 0.39 +Other | | 0.067 | | | 0.08 + +Nlocal: 900 ave 900 max 900 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 762 ave 762 max 762 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 40697 ave 40697 max 40697 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 40697 +Ave neighs/atom = 45.2189 +Neighbor list builds = 32 +Dangerous builds = 0 +Total wall time: 0:01:20 diff --git a/examples/USER/sdpd/2d-diffusion/log.24Oct18.2d-diffusion.g++.4 b/examples/USER/sdpd/2d-diffusion/log.24Oct18.2d-diffusion.g++.4 new file mode 100644 index 0000000000..f904b78ab4 --- /dev/null +++ b/examples/USER/sdpd/2d-diffusion/log.24Oct18.2d-diffusion.g++.4 @@ -0,0 +1,226 @@ +LAMMPS (24 Oct 2018) +dimension 2 +units micro +atom_style meso + +variable R equal 0.5 # radius of sphere micrometers +variable a equal $R/5 # lattice spacing micrometers +variable a equal 0.5/5 +variable L equal $R*3 +variable L equal 0.5*3 +variable T equal 300. +variable rho_0 equal 1. # density picograms/micrometer^3 +variable c_0 equal 100. # speed of sound micrometers/microsecond +variable mu equal 1. # dynamic viscosity picogram/(micrometer-microsecond) +variable h equal $a*4.5 # kernel function cutoff micrometers +variable h equal 0.1*4.5 +variable mass equal $a*$a*$a*${rho_0} +variable mass equal 0.1*$a*$a*${rho_0} +variable mass equal 0.1*0.1*$a*${rho_0} +variable mass equal 0.1*0.1*0.1*${rho_0} +variable mass equal 0.1*0.1*0.1*1 +variable dt equal 1e-3 # timestep microseconds +variable skin equal 0.2*$h +variable skin equal 0.2*0.45 + +region box block -$L $L -$L $L 0 $a units box +region box block -1.5 $L -$L $L 0 $a units box +region box block -1.5 1.5 -$L $L 0 $a units box +region box block -1.5 1.5 -1.5 $L 0 $a units box +region box block -1.5 1.5 -1.5 1.5 0 $a units box +region box block -1.5 1.5 -1.5 1.5 0 0.1 units box +create_box 2 box +Created orthogonal box = (-1.5 -1.5 0) to (1.5 1.5 0.1) + 2 by 2 by 1 MPI processor grid +lattice sq $a +lattice sq 0.1 +Lattice spacing in x,y,z = 0.1 0.1 0.1 + +create_atoms 1 box +Created 900 atoms + Time spent = 0.0010246 secs + +region sphere sphere 0 0 0 $R units box +region sphere sphere 0 0 0 0.5 units box +set region sphere type 2 + 81 settings made for type + +group fluid type 1 +819 atoms in group fluid +group sphere type 2 +81 atoms in group sphere + +mass * ${mass} +mass * 0.001 +set group all meso/rho ${rho_0} +set group all meso/rho 1 + 900 settings made for meso/rho + +pair_style sdpd/taitwater/isothermal $T ${mu} 76787 # temperature viscosity random_seed +pair_style sdpd/taitwater/isothermal 300 ${mu} 76787 +pair_style sdpd/taitwater/isothermal 300 1 76787 +pair_coeff * * ${rho_0} ${c_0} ${h} +pair_coeff * * 1 ${c_0} ${h} +pair_coeff * * 1 100 ${h} +pair_coeff * * 1 100 0.45 + +fix 1 fluid meso +fix 2 sphere rigid/meso single +1 rigid bodies with 81 atoms + +fix 2d all enforce2d + +neighbor ${skin} bin +neighbor 0.09 bin +neigh_modify delay 0 every 1 check yes +timestep ${dt} +timestep 0.001 + +dump dump_id all atom 100 dump.lammpstrj + +thermo 100 +thermo_style custom step time nbuild ndanger + +run 10000 +Neighbor list info ... + update every 1 steps, delay 0 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 0.54 + ghost atom cutoff = 0.54 + binsize = 0.27, bins = 12 12 1 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair sdpd/taitwater/isothermal, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/2d/newton + bin: standard +Per MPI rank memory allocation (min/avg/max) = 6.087 | 6.087 | 6.087 Mbytes +Step Time Nbuild Ndanger + 0 0 0 0 + 100 0.1 0 0 + 200 0.2 0 0 + 300 0.3 0 0 + 400 0.4 1 0 + 500 0.5 1 0 + 600 0.6 1 0 + 700 0.7 2 0 + 800 0.8 2 0 + 900 0.9 2 0 + 1000 1 3 0 + 1100 1.1 3 0 + 1200 1.2 3 0 + 1300 1.3 4 0 + 1400 1.4 4 0 + 1500 1.5 5 0 + 1600 1.6 5 0 + 1700 1.7 5 0 + 1800 1.8 6 0 + 1900 1.9 6 0 + 2000 2 6 0 + 2100 2.1 7 0 + 2200 2.2 7 0 + 2300 2.3 7 0 + 2400 2.4 8 0 + 2500 2.5 8 0 + 2600 2.6 8 0 + 2700 2.7 9 0 + 2800 2.8 9 0 + 2900 2.9 9 0 + 3000 3 9 0 + 3100 3.1 10 0 + 3200 3.2 10 0 + 3300 3.3 10 0 + 3400 3.4 11 0 + 3500 3.5 11 0 + 3600 3.6 11 0 + 3700 3.7 12 0 + 3800 3.8 12 0 + 3900 3.9 12 0 + 4000 4 13 0 + 4100 4.1 13 0 + 4200 4.2 13 0 + 4300 4.3 14 0 + 4400 4.4 14 0 + 4500 4.5 15 0 + 4600 4.6 15 0 + 4700 4.7 15 0 + 4800 4.8 16 0 + 4900 4.9 16 0 + 5000 5 16 0 + 5100 5.1 16 0 + 5200 5.2 17 0 + 5300 5.3 17 0 + 5400 5.4 18 0 + 5500 5.5 18 0 + 5600 5.6 19 0 + 5700 5.7 19 0 + 5800 5.8 19 0 + 5900 5.9 20 0 + 6000 6 20 0 + 6100 6.1 20 0 + 6200 6.2 21 0 + 6300 6.3 21 0 + 6400 6.4 21 0 + 6500 6.5 22 0 + 6600 6.6 22 0 + 6700 6.7 22 0 + 6800 6.8 23 0 + 6900 6.9 23 0 + 7000 7 23 0 + 7100 7.1 24 0 + 7200 7.2 24 0 + 7300 7.3 24 0 + 7400 7.4 25 0 + 7500 7.5 25 0 + 7600 7.6 25 0 + 7700 7.7 26 0 + 7800 7.8 26 0 + 7900 7.9 26 0 + 8000 8 27 0 + 8100 8.1 27 0 + 8200 8.2 27 0 + 8300 8.3 28 0 + 8400 8.4 28 0 + 8500 8.5 28 0 + 8600 8.6 28 0 + 8700 8.7 29 0 + 8800 8.8 29 0 + 8900 8.9 29 0 + 9000 9 30 0 + 9100 9.1 30 0 + 9200 9.2 31 0 + 9300 9.3 31 0 + 9400 9.4 31 0 + 9500 9.5 31 0 + 9600 9.6 32 0 + 9700 9.7 32 0 + 9800 9.8 32 0 + 9900 9.9 33 0 + 10000 10 33 0 +Loop time of 69.01 on 4 procs for 10000 steps with 900 atoms + +Performance: 12519931.275 ns/day, 0.000 hours/ns, 144.907 timesteps/s +48.7% CPU use with 4 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 56.528 | 57.936 | 58.729 | 11.0 | 83.95 +Neigh | 0.013157 | 0.013382 | 0.013551 | 0.1 | 0.02 +Comm | 8.9594 | 9.7555 | 11.113 | 26.7 | 14.14 +Output | 0.14644 | 0.15009 | 0.15809 | 1.2 | 0.22 +Modify | 0.72913 | 0.91574 | 1.0524 | 12.4 | 1.33 +Other | | 0.2389 | | | 0.35 + +Nlocal: 225 ave 229 max 223 min +Histogram: 1 2 0 0 0 0 0 0 0 1 +Nghost: 442 ave 444 max 439 min +Histogram: 1 0 0 0 1 0 0 0 0 2 +Neighs: 10188.8 ave 10437 max 9932 min +Histogram: 1 0 0 1 0 0 0 1 0 1 + +Total # of neighbors = 40755 +Ave neighs/atom = 45.2833 +Neighbor list builds = 33 +Dangerous builds = 0 +Total wall time: 0:01:09 diff --git a/examples/USER/sdpd/README b/examples/USER/sdpd/README new file mode 100644 index 0000000000..72feff15dc --- /dev/null +++ b/examples/USER/sdpd/README @@ -0,0 +1,24 @@ +Smoothed Dissipative Particle Dynamics examples + +equipartition-verification: + This example verifies the equipartition theorem. + It simulates a periodic box of water with no solid bodies. + If equipartition theorem holds true, the average of each component of + translational kinetic energy should be equal to k_B T, and therefore + vx_sq_check, vy_sq_check, and vz_sq_check should fluctuate near 1. + +2d-diffusion: + This example demonstrates the free diffusion of a disk in 2D. + The 3D simulation is similar but takes much longer to complete. + As with other statistical experiments you need an ensemble to + extract meaningful average quantities. + For a more realistic simulation you should increase the resolution + of the disk/sphere which also necessitates reduction of timestep. + +2d-diffusion-in-shear-flow: + This example demonstrates the diffusion of a disk in shear flow in 2D. + The 3D simulation is similar but takes much longer to complete. + As with other statistical experiments you need an ensemble to + extract meaningful average quantities. + For a more realistic simulation you should increase the resolution + of the disk/sphere which also necessitates reduction of timestep. diff --git a/examples/USER/sdpd/equipartition-verification/in.lammps b/examples/USER/sdpd/equipartition-verification/in.lammps new file mode 100644 index 0000000000..0d06723f59 --- /dev/null +++ b/examples/USER/sdpd/equipartition-verification/in.lammps @@ -0,0 +1,45 @@ +dimension 3 +units micro +atom_style meso + +variable a equal 0.1 # lattice spacing micrometers +variable L equal $a*10 +variable T equal 300. +variable kB equal 1.3806504e-8 # picogram-micrometer^2/(microsecond^2-Kelvin) +variable rho_0 equal 1. # density picograms/micrometer^3 +variable c_0 equal 10. # speed of sound micrometers/microsecond +variable mu equal 1. # dynamic viscosity picogram/(micrometer-microsecond) +variable h equal $a*4.0 # kernel function cutoff micrometers +variable mass equal $a*$a*$a*${rho_0} +variable dt equal 5e-4 # timestep microseconds +variable skin equal 0.1*$h + +region box block -$L $L -$L $L -$L $L units box +create_box 1 box +lattice sc $a + +create_atoms 1 box + +mass * ${mass} +set group all meso/rho ${rho_0} + +pair_style sdpd/taitwater/isothermal $T ${mu} 76787 # temperature viscosity random_seed +pair_coeff * * ${rho_0} ${c_0} ${h} + +variable vx_sq atom vx*vx +variable vy_sq atom vy*vy +variable vz_sq atom vz*vz +compute v_sq all reduce ave v_vx_sq v_vy_sq v_vz_sq +variable vx_sq_check equal c_v_sq[1]*${mass}/${kB}/$T +variable vy_sq_check equal c_v_sq[2]*${mass}/${kB}/$T +variable vz_sq_check equal c_v_sq[3]*${mass}/${kB}/$T + +fix 1 all meso + +neighbor ${skin} bin +timestep ${dt} + +thermo 10 +thermo_style custom step time v_vx_sq_check v_vy_sq_check v_vz_sq_check + +run 200 diff --git a/examples/USER/sdpd/equipartition-verification/log.24Oct18.equipartition.g++.1 b/examples/USER/sdpd/equipartition-verification/log.24Oct18.equipartition.g++.1 new file mode 100644 index 0000000000..06ffd699bc --- /dev/null +++ b/examples/USER/sdpd/equipartition-verification/log.24Oct18.equipartition.g++.1 @@ -0,0 +1,146 @@ +LAMMPS (24 Oct 2018) +dimension 3 +units micro +atom_style meso + +variable a equal 0.1 # lattice spacing micrometers +variable L equal $a*10 +variable L equal 0.1*10 +variable T equal 300. +variable kB equal 1.3806504e-8 # picogram-micrometer^2/(microsecond^2-Kelvin) +variable rho_0 equal 1. # density picograms/micrometer^3 +variable c_0 equal 10. # speed of sound micrometers/microsecond +variable mu equal 1. # dynamic viscosity picogram/(micrometer-microsecond) +variable h equal $a*4.0 # kernel function cutoff micrometers +variable h equal 0.1*4.0 +variable mass equal $a*$a*$a*${rho_0} +variable mass equal 0.1*$a*$a*${rho_0} +variable mass equal 0.1*0.1*$a*${rho_0} +variable mass equal 0.1*0.1*0.1*${rho_0} +variable mass equal 0.1*0.1*0.1*1 +variable dt equal 5e-4 # timestep microseconds +variable skin equal 0.1*$h +variable skin equal 0.1*0.4 + +region box block -$L $L -$L $L -$L $L units box +region box block -1 $L -$L $L -$L $L units box +region box block -1 1 -$L $L -$L $L units box +region box block -1 1 -1 $L -$L $L units box +region box block -1 1 -1 1 -$L $L units box +region box block -1 1 -1 1 -1 $L units box +region box block -1 1 -1 1 -1 1 units box +create_box 1 box +Created orthogonal box = (-1 -1 -1) to (1 1 1) + 1 by 1 by 1 MPI processor grid +lattice sc $a +lattice sc 0.1 +Lattice spacing in x,y,z = 0.1 0.1 0.1 + +create_atoms 1 box +Created 8000 atoms + Time spent = 0.00285411 secs + +mass * ${mass} +mass * 0.001 +set group all meso/rho ${rho_0} +set group all meso/rho 1 + 8000 settings made for meso/rho + +pair_style sdpd/taitwater/isothermal $T ${mu} 76787 # temperature viscosity random_seed +pair_style sdpd/taitwater/isothermal 300 ${mu} 76787 +pair_style sdpd/taitwater/isothermal 300 1 76787 +pair_coeff * * ${rho_0} ${c_0} ${h} +pair_coeff * * 1 ${c_0} ${h} +pair_coeff * * 1 10 ${h} +pair_coeff * * 1 10 0.4 + +variable vx_sq atom vx*vx +variable vy_sq atom vy*vy +variable vz_sq atom vz*vz +compute v_sq all reduce ave v_vx_sq v_vy_sq v_vz_sq +variable vx_sq_check equal c_v_sq[1]*${mass}/${kB}/$T +variable vx_sq_check equal c_v_sq[1]*0.001/${kB}/$T +variable vx_sq_check equal c_v_sq[1]*0.001/1.3806504e-08/$T +variable vx_sq_check equal c_v_sq[1]*0.001/1.3806504e-08/300 +variable vy_sq_check equal c_v_sq[2]*${mass}/${kB}/$T +variable vy_sq_check equal c_v_sq[2]*0.001/${kB}/$T +variable vy_sq_check equal c_v_sq[2]*0.001/1.3806504e-08/$T +variable vy_sq_check equal c_v_sq[2]*0.001/1.3806504e-08/300 +variable vz_sq_check equal c_v_sq[3]*${mass}/${kB}/$T +variable vz_sq_check equal c_v_sq[3]*0.001/${kB}/$T +variable vz_sq_check equal c_v_sq[3]*0.001/1.3806504e-08/$T +variable vz_sq_check equal c_v_sq[3]*0.001/1.3806504e-08/300 + +fix 1 all meso + +neighbor ${skin} bin +neighbor 0.04 bin +timestep ${dt} +timestep 0.0005 + +thermo 10 +thermo_style custom step time v_vx_sq_check v_vy_sq_check v_vz_sq_check + +run 200 +Neighbor list info ... + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 0.44 + ghost atom cutoff = 0.44 + binsize = 0.22, bins = 10 10 10 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair sdpd/taitwater/isothermal, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/3d/newton + bin: standard +Per MPI rank memory allocation (min/avg/max) = 13.54 | 13.54 | 13.54 Mbytes +Step Time v_vx_sq_check v_vy_sq_check v_vz_sq_check + 0 0 0 0 0 + 10 0.005 0.70973271 0.71495693 0.71910087 + 20 0.01 0.90418096 0.88845437 0.89659567 + 30 0.015 0.9590736 0.97880338 0.9619016 + 40 0.02 0.98533774 0.96057682 0.95600448 + 50 0.025 0.96433662 0.96650071 0.95509683 + 60 0.03 0.96598029 0.96373656 0.96734888 + 70 0.035 0.95433045 0.98004764 0.96255924 + 80 0.04 0.97872906 0.95987289 0.96623598 + 90 0.045 0.99913888 0.99255731 0.95616142 + 100 0.05 0.98872675 0.97141018 0.95338841 + 110 0.055 0.97794592 0.97389258 0.98473719 + 120 0.06 0.98389266 0.96716284 0.95504862 + 130 0.065 0.98572886 0.96680923 0.95599065 + 140 0.07 0.97602684 0.97580081 0.9886878 + 150 0.075 0.99172003 0.95027467 0.96028033 + 160 0.08 0.96793247 0.94590928 0.95644301 + 170 0.085 0.94167619 0.98048861 0.93439426 + 180 0.09 0.97277934 0.97383622 0.96900866 + 190 0.095 0.96647288 1.0027643 0.96230782 + 200 0.1 0.94864291 0.95902585 0.96398175 +Loop time of 60.1095 on 1 procs for 200 steps with 8000 atoms + +Performance: 143737.595 ns/day, 0.000 hours/ns, 3.327 timesteps/s +99.7% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 59.92 | 59.92 | 59.92 | 0.0 | 99.68 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 0.11154 | 0.11154 | 0.11154 | 0.0 | 0.19 +Output | 0.0063498 | 0.0063498 | 0.0063498 | 0.0 | 0.01 +Modify | 0.043546 | 0.043546 | 0.043546 | 0.0 | 0.07 +Other | | 0.02811 | | | 0.05 + +Nlocal: 8000 ave 8000 max 8000 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 16389 ave 16389 max 16389 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 1.456e+06 ave 1.456e+06 max 1.456e+06 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 1456000 +Ave neighs/atom = 182 +Neighbor list builds = 0 +Dangerous builds = 0 +Total wall time: 0:01:00 diff --git a/examples/USER/sdpd/equipartition-verification/log.24Oct18.equipartition.g++.4 b/examples/USER/sdpd/equipartition-verification/log.24Oct18.equipartition.g++.4 new file mode 100644 index 0000000000..88509f0fd1 --- /dev/null +++ b/examples/USER/sdpd/equipartition-verification/log.24Oct18.equipartition.g++.4 @@ -0,0 +1,146 @@ +LAMMPS (24 Oct 2018) +dimension 3 +units micro +atom_style meso + +variable a equal 0.1 # lattice spacing micrometers +variable L equal $a*10 +variable L equal 0.1*10 +variable T equal 300. +variable kB equal 1.3806504e-8 # picogram-micrometer^2/(microsecond^2-Kelvin) +variable rho_0 equal 1. # density picograms/micrometer^3 +variable c_0 equal 10. # speed of sound micrometers/microsecond +variable mu equal 1. # dynamic viscosity picogram/(micrometer-microsecond) +variable h equal $a*4.0 # kernel function cutoff micrometers +variable h equal 0.1*4.0 +variable mass equal $a*$a*$a*${rho_0} +variable mass equal 0.1*$a*$a*${rho_0} +variable mass equal 0.1*0.1*$a*${rho_0} +variable mass equal 0.1*0.1*0.1*${rho_0} +variable mass equal 0.1*0.1*0.1*1 +variable dt equal 5e-4 # timestep microseconds +variable skin equal 0.1*$h +variable skin equal 0.1*0.4 + +region box block -$L $L -$L $L -$L $L units box +region box block -1 $L -$L $L -$L $L units box +region box block -1 1 -$L $L -$L $L units box +region box block -1 1 -1 $L -$L $L units box +region box block -1 1 -1 1 -$L $L units box +region box block -1 1 -1 1 -1 $L units box +region box block -1 1 -1 1 -1 1 units box +create_box 1 box +Created orthogonal box = (-1 -1 -1) to (1 1 1) + 1 by 2 by 2 MPI processor grid +lattice sc $a +lattice sc 0.1 +Lattice spacing in x,y,z = 0.1 0.1 0.1 + +create_atoms 1 box +Created 8000 atoms + Time spent = 0.00252754 secs + +mass * ${mass} +mass * 0.001 +set group all meso/rho ${rho_0} +set group all meso/rho 1 + 8000 settings made for meso/rho + +pair_style sdpd/taitwater/isothermal $T ${mu} 76787 # temperature viscosity random_seed +pair_style sdpd/taitwater/isothermal 300 ${mu} 76787 +pair_style sdpd/taitwater/isothermal 300 1 76787 +pair_coeff * * ${rho_0} ${c_0} ${h} +pair_coeff * * 1 ${c_0} ${h} +pair_coeff * * 1 10 ${h} +pair_coeff * * 1 10 0.4 + +variable vx_sq atom vx*vx +variable vy_sq atom vy*vy +variable vz_sq atom vz*vz +compute v_sq all reduce ave v_vx_sq v_vy_sq v_vz_sq +variable vx_sq_check equal c_v_sq[1]*${mass}/${kB}/$T +variable vx_sq_check equal c_v_sq[1]*0.001/${kB}/$T +variable vx_sq_check equal c_v_sq[1]*0.001/1.3806504e-08/$T +variable vx_sq_check equal c_v_sq[1]*0.001/1.3806504e-08/300 +variable vy_sq_check equal c_v_sq[2]*${mass}/${kB}/$T +variable vy_sq_check equal c_v_sq[2]*0.001/${kB}/$T +variable vy_sq_check equal c_v_sq[2]*0.001/1.3806504e-08/$T +variable vy_sq_check equal c_v_sq[2]*0.001/1.3806504e-08/300 +variable vz_sq_check equal c_v_sq[3]*${mass}/${kB}/$T +variable vz_sq_check equal c_v_sq[3]*0.001/${kB}/$T +variable vz_sq_check equal c_v_sq[3]*0.001/1.3806504e-08/$T +variable vz_sq_check equal c_v_sq[3]*0.001/1.3806504e-08/300 + +fix 1 all meso + +neighbor ${skin} bin +neighbor 0.04 bin +timestep ${dt} +timestep 0.0005 + +thermo 10 +thermo_style custom step time v_vx_sq_check v_vy_sq_check v_vz_sq_check + +run 200 +Neighbor list info ... + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 0.44 + ghost atom cutoff = 0.44 + binsize = 0.22, bins = 10 10 10 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair sdpd/taitwater/isothermal, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/3d/newton + bin: standard +Per MPI rank memory allocation (min/avg/max) = 5.795 | 5.795 | 5.795 Mbytes +Step Time v_vx_sq_check v_vy_sq_check v_vz_sq_check + 0 0 0 0 0 + 10 0.005 0.71224819 0.71470372 0.7008956 + 20 0.01 0.90627589 0.90683966 0.90116506 + 30 0.015 0.938505 0.95884272 0.93337542 + 40 0.02 0.94394649 0.93668038 0.96468004 + 50 0.025 0.97152309 0.97546161 0.95107762 + 60 0.03 0.94710871 0.95678322 0.97285504 + 70 0.035 0.96253148 0.95838642 0.95450883 + 80 0.04 0.97581495 0.95278681 0.95099478 + 90 0.045 0.96251614 0.9740684 0.96081505 + 100 0.05 0.94191275 0.97137523 0.94084858 + 110 0.055 0.953406 0.95739684 0.98574522 + 120 0.06 0.99001614 0.99608287 0.9839996 + 130 0.065 0.96575225 0.94309655 0.92847798 + 140 0.07 0.97642687 0.97458638 0.94696406 + 150 0.075 0.99316381 0.96876814 0.95440106 + 160 0.08 0.94589744 0.95264791 0.95495169 + 170 0.085 0.97599092 0.95336014 0.97687718 + 180 0.09 0.97214242 0.9726305 0.9726035 + 190 0.095 0.97577583 0.96523645 0.9756968 + 200 0.1 0.96386053 0.97268854 0.94582436 +Loop time of 32.5247 on 4 procs for 200 steps with 8000 atoms + +Performance: 265644.515 ns/day, 0.000 hours/ns, 6.149 timesteps/s +73.9% CPU use with 4 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 27.385 | 28.409 | 28.761 | 11.1 | 87.34 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 3.582 | 3.9343 | 4.9531 | 29.7 | 12.10 +Output | 0.022267 | 0.026073 | 0.033141 | 2.7 | 0.08 +Modify | 0.031714 | 0.033134 | 0.034367 | 0.6 | 0.10 +Other | | 0.1226 | | | 0.38 + +Nlocal: 2000 ave 2000 max 2000 min +Histogram: 4 0 0 0 0 0 0 0 0 0 +Nghost: 8469 ave 8469 max 8469 min +Histogram: 4 0 0 0 0 0 0 0 0 0 +Neighs: 364000 ave 376628 max 351184 min +Histogram: 1 0 1 0 0 0 0 1 0 1 + +Total # of neighbors = 1456000 +Ave neighs/atom = 182 +Neighbor list builds = 0 +Dangerous builds = 0 +Total wall time: 0:00:32 diff --git a/src/Depend.sh b/src/Depend.sh index a3f4e49667..b7ee740c52 100755 --- a/src/Depend.sh +++ b/src/Depend.sh @@ -107,6 +107,7 @@ fi if (test $1 = "RIGID") then depend USER-OMP + depend USER-SDPD fi if (test $1 = "SNAP") then diff --git a/src/USER-SDPD/Install.sh b/src/USER-SDPD/Install.sh new file mode 100644 index 0000000000..cc4af1c038 --- /dev/null +++ b/src/USER-SDPD/Install.sh @@ -0,0 +1,37 @@ +# Install/Uninstall package files in LAMMPS +# mode = 0/1/2 for uninstall/install/update + +mode=$1 + +# enforce using portable C locale +LC_ALL=C +export LC_ALL + +# arg1 = file, arg2 = file it depends on + +action () { + if (test $mode = 0) then + rm -f ../$1 + elif (! cmp -s $1 ../$1) then + if (test -z "$2" || test -e ../$2) then + cp $1 .. + if (test $mode = 2) then + echo " updating src/$1" + fi + fi + elif (test -n "$2") then + if (test ! -e ../$2) then + rm -f ../$1 + fi + fi +} + +# package files without dependencies +action pair_sdpd_taitwater_isothermal.h +action pair_sdpd_taitwater_isothermal.cpp +action fix_meso_move.h +action fix_meso_move.cpp + +# package files with dependencies +action fix_rigid_meso.h fix_rigid.h +action fix_rigid_meso.cpp fix_rigid.h diff --git a/src/USER-SDPD/README b/src/USER-SDPD/README new file mode 100644 index 0000000000..88ece7ddbf --- /dev/null +++ b/src/USER-SDPD/README @@ -0,0 +1,13 @@ +This package implements the Smoothed Dissipative Particle Dynamics (SDPD) +method for modelling fluids at mesoscale and diffusion of colloids. +Currently it adds pair style sdpd/taitwater/isothermal for modelling water +at isothermal conditions using the Tait equation of state. +It also adds fix meso/move command to move mesoscopic SPH/SDPD particles with +prescribed velocity and fix rigid/meso command to integrate rigid bodies +composed of mesoscopic SPH/SDPD particles. + +Creator of this package is: + +Morteza Jalalvand +Institute for Advanced Studies in Basic Sciences +jalalvand.m AT gmail DOT com diff --git a/src/USER-SDPD/fix_meso_move.cpp b/src/USER-SDPD/fix_meso_move.cpp new file mode 100644 index 0000000000..3e2902d0cc --- /dev/null +++ b/src/USER-SDPD/fix_meso_move.cpp @@ -0,0 +1,998 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: + Morteza Jalalvand (IASBS) jalalvand.m AT gmail.com +------------------------------------------------------------------------- */ + +#include "string.h" +#include "math.h" +#include "fix_meso_move.h" +#include "atom.h" +#include "group.h" +#include "update.h" +#include "modify.h" +#include "force.h" +#include "domain.h" +#include "lattice.h" +#include "comm.h" +#include "input.h" +#include "variable.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace MathConst; + +enum{LINEAR,WIGGLE,ROTATE,VARIABLE}; +enum{EQUAL,ATOM}; + +/* ---------------------------------------------------------------------- */ + +FixMesoMove::FixMesoMove (LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), + xvarstr(NULL), yvarstr(NULL), zvarstr(NULL), + vxvarstr(NULL), vyvarstr(NULL), vzvarstr(NULL), + xoriginal(NULL), displace(NULL), velocity(NULL) { + if ((atom->e_flag != 1) || (atom->rho_flag != 1)) + error->all(FLERR, + "fix meso/move command requires atom_style with both energy and density"); + + if (narg < 4) error->all(FLERR,"Illegal fix meso/move command"); + + restart_global = 1; + restart_peratom = 1; + peratom_flag = 1; + size_peratom_cols = 3; + peratom_freq = 1; + time_integrate = 1; + create_attribute = 1; + displaceflag = 0; + velocityflag = 0; + maxatom = 0; + + // parse args + + int iarg; + + if (strcmp(arg[3],"linear") == 0) { + if (narg < 7) error->all(FLERR,"Illegal fix meso/move command"); + iarg = 7; + mstyle = LINEAR; + if (strcmp(arg[4],"NULL") == 0) vxflag = 0; + else { + vxflag = 1; + vx = force->numeric(FLERR,arg[4]); + } + if (strcmp(arg[5],"NULL") == 0) vyflag = 0; + else { + vyflag = 1; + vy = force->numeric(FLERR,arg[5]); + } + if (strcmp(arg[6],"NULL") == 0) vzflag = 0; + else { + vzflag = 1; + vz = force->numeric(FLERR,arg[6]); + } + + } else if (strcmp(arg[3],"wiggle") == 0) { + if (narg < 8) error->all(FLERR,"Illegal fix meso/move command"); + iarg = 8; + mstyle = WIGGLE; + if (strcmp(arg[4],"NULL") == 0) axflag = 0; + else { + axflag = 1; + ax = force->numeric(FLERR,arg[4]); + } + if (strcmp(arg[5],"NULL") == 0) ayflag = 0; + else { + ayflag = 1; + ay = force->numeric(FLERR,arg[5]); + } + if (strcmp(arg[6],"NULL") == 0) azflag = 0; + else { + azflag = 1; + az = force->numeric(FLERR,arg[6]); + } + period = force->numeric(FLERR,arg[7]); + if (period <= 0.0) error->all(FLERR,"Illegal fix meso/move command"); + + } else if (strcmp(arg[3],"rotate") == 0) { + if (narg < 11) error->all(FLERR,"Illegal fix meso/move command"); + iarg = 11; + mstyle = ROTATE; + point[0] = force->numeric(FLERR,arg[4]); + point[1] = force->numeric(FLERR,arg[5]); + point[2] = force->numeric(FLERR,arg[6]); + axis[0] = force->numeric(FLERR,arg[7]); + axis[1] = force->numeric(FLERR,arg[8]); + axis[2] = force->numeric(FLERR,arg[9]); + period = force->numeric(FLERR,arg[10]); + if (period <= 0.0) error->all(FLERR,"Illegal fix meso/move command"); + + } else if (strcmp(arg[3],"variable") == 0) { + if (narg < 10) error->all(FLERR,"Illegal fix meso/move command"); + iarg = 10; + mstyle = VARIABLE; + if (strcmp(arg[4],"NULL") == 0) xvarstr = NULL; + else if (strstr(arg[4],"v_") == arg[4]) { + int n = strlen(&arg[4][2]) + 1; + xvarstr = new char[n]; + strcpy(xvarstr,&arg[4][2]); + } else error->all(FLERR,"Illegal fix meso/move command"); + if (strcmp(arg[5],"NULL") == 0) yvarstr = NULL; + else if (strstr(arg[5],"v_") == arg[5]) { + int n = strlen(&arg[5][2]) + 1; + yvarstr = new char[n]; + strcpy(yvarstr,&arg[5][2]); + } else error->all(FLERR,"Illegal fix meso/move command"); + if (strcmp(arg[6],"NULL") == 0) zvarstr = NULL; + else if (strstr(arg[6],"v_") == arg[6]) { + int n = strlen(&arg[6][2]) + 1; + zvarstr = new char[n]; + strcpy(zvarstr,&arg[6][2]); + } else error->all(FLERR,"Illegal fix meso/move command"); + if (strcmp(arg[7],"NULL") == 0) vxvarstr = NULL; + else if (strstr(arg[7],"v_") == arg[7]) { + int n = strlen(&arg[7][2]) + 1; + vxvarstr = new char[n]; + strcpy(vxvarstr,&arg[7][2]); + } else error->all(FLERR,"Illegal fix meso/move command"); + if (strcmp(arg[8],"NULL") == 0) vyvarstr = NULL; + else if (strstr(arg[8],"v_") == arg[8]) { + int n = strlen(&arg[8][2]) + 1; + vyvarstr = new char[n]; + strcpy(vyvarstr,&arg[8][2]); + } else error->all(FLERR,"Illegal fix meso/move command"); + if (strcmp(arg[9],"NULL") == 0) vzvarstr = NULL; + else if (strstr(arg[9],"v_") == arg[9]) { + int n = strlen(&arg[9][2]) + 1; + vzvarstr = new char[n]; + strcpy(vzvarstr,&arg[9][2]); + } else error->all(FLERR,"Illegal fix meso/move command"); + + } else error->all(FLERR,"Illegal fix meso/move command"); + + // optional args + + int scaleflag = 1; + + while (iarg < narg) { + if (strcmp(arg[iarg],"units") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix meso/move command"); + if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; + else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; + else error->all(FLERR,"Illegal fix meso/move command"); + iarg += 2; + } else error->all(FLERR,"Illegal fix meso/move command"); + } + + // error checks and warnings + + if (domain->dimension == 2) { + if (mstyle == LINEAR && vzflag && vz != 0.0) + error->all(FLERR,"Fix meso/move cannot set linear z motion for 2d problem"); + if (mstyle == WIGGLE && azflag && az != 0.0) + error->all(FLERR,"Fix meso/move cannot set wiggle z motion for 2d problem"); + if (mstyle == ROTATE && (axis[0] != 0.0 || axis[1] != 0.0)) + error->all(FLERR, "Fix meso/move cannot rotate aroung non z-axis for 2d problem"); + if (mstyle == VARIABLE && (zvarstr || vzvarstr)) + error->all(FLERR, "Fix meso/move cannot define z or vz variable for 2d problem"); + } + + // setup scaling and apply scaling factors to velocity & amplitude + + if ((mstyle == LINEAR || mstyle == WIGGLE || mstyle == ROTATE) && + scaleflag) { + double xscale = domain->lattice->xlattice; + double yscale = domain->lattice->ylattice; + double zscale = domain->lattice->zlattice; + + if (mstyle == LINEAR) { + if (vxflag) vx *= xscale; + if (vyflag) vy *= yscale; + if (vzflag) vz *= zscale; + } else if (mstyle == WIGGLE) { + if (axflag) ax *= xscale; + if (ayflag) ay *= yscale; + if (azflag) az *= zscale; + } else if (mstyle == ROTATE) { + point[0] *= xscale; + point[1] *= yscale; + point[2] *= zscale; + } + } + + // set omega_rotate from period + + if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = MY_2PI / period; + + // runit = unit vector along rotation axis + + if (mstyle == ROTATE) { + double len = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]); + if (len == 0.0) + error->all(FLERR,"Zero length rotation vector with fix meso/move"); + runit[0] = axis[0]/len; + runit[1] = axis[1]/len; + runit[2] = axis[2]/len; + } + + // perform initial allocation of atom-based array + // register with Atom class + + grow_arrays(atom->nmax); + atom->add_callback(0); + atom->add_callback(1); + + displace = velocity = NULL; + + // xoriginal = initial unwrapped positions of atoms + + double **x = atom->x; + imageint *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]); + else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; + } + + time_origin = update->ntimestep; +} + +/* ---------------------------------------------------------------------- */ + +FixMesoMove::~FixMesoMove () { + // unregister callbacks to this fix from Atom class + + atom->delete_callback(id,0); + atom->delete_callback(id,1); + + // delete locally stored arrays + + memory->destroy(xoriginal); + memory->destroy(displace); + memory->destroy(velocity); + + delete [] xvarstr; + delete [] yvarstr; + delete [] zvarstr; + delete [] vxvarstr; + delete [] vyvarstr; + delete [] vzvarstr; +} + +/* ---------------------------------------------------------------------- */ + +int FixMesoMove::setmask () { + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + mask |= PRE_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixMesoMove::init () { + dt = update->dt; + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + + // set indices and style of all variables + + displaceflag = velocityflag = 0; + + if (mstyle == VARIABLE) { + if (xvarstr) { + xvar = input->variable->find(xvarstr); + if (xvar < 0) error->all(FLERR, "Variable name for fix meso/move does not exist"); + if (input->variable->equalstyle(xvar)) xvarstyle = EQUAL; + else if (input->variable->atomstyle(xvar)) xvarstyle = ATOM; + else error->all(FLERR,"Variable for fix meso/move is invalid style"); + } + if (yvarstr) { + yvar = input->variable->find(yvarstr); + if (yvar < 0) error->all(FLERR, "Variable name for fix meso/move does not exist"); + if (input->variable->equalstyle(yvar)) yvarstyle = EQUAL; + else if (input->variable->atomstyle(yvar)) yvarstyle = ATOM; + else error->all(FLERR,"Variable for fix meso/move is invalid style"); + } + if (zvarstr) { + zvar = input->variable->find(zvarstr); + if (zvar < 0) error->all(FLERR, "Variable name for fix meso/move does not exist"); + if (input->variable->equalstyle(zvar)) zvarstyle = EQUAL; + else if (input->variable->atomstyle(zvar)) zvarstyle = ATOM; + else error->all(FLERR,"Variable for fix meso/move is invalid style"); + } + if (vxvarstr) { + vxvar = input->variable->find(vxvarstr); + if (vxvar < 0) error->all(FLERR, "Variable name for fix meso/move does not exist"); + if (input->variable->equalstyle(vxvar)) vxvarstyle = EQUAL; + else if (input->variable->atomstyle(vxvar)) vxvarstyle = ATOM; + else error->all(FLERR,"Variable for fix meso/move is invalid style"); + } + if (vyvarstr) { + vyvar = input->variable->find(vyvarstr); + if (vyvar < 0) error->all(FLERR, "Variable name for fix meso/move does not exist"); + if (input->variable->equalstyle(vyvar)) vyvarstyle = EQUAL; + else if (input->variable->atomstyle(vyvar)) vyvarstyle = ATOM; + else error->all(FLERR,"Variable for fix meso/move is invalid style"); + } + if (vzvarstr) { + vzvar = input->variable->find(vzvarstr); + if (vzvar < 0) error->all(FLERR, "Variable name for fix meso/move does not exist"); + if (input->variable->equalstyle(vzvar)) vzvarstyle = EQUAL; + else if (input->variable->atomstyle(vzvar)) vzvarstyle = ATOM; + else error->all(FLERR,"Variable for fix meso/move is invalid style"); + } + + if (xvarstr && xvarstyle == ATOM) displaceflag = 1; + if (yvarstr && yvarstyle == ATOM) displaceflag = 1; + if (zvarstr && zvarstyle == ATOM) displaceflag = 1; + if (vxvarstr && vxvarstyle == ATOM) velocityflag = 1; + if (vyvarstr && vyvarstyle == ATOM) velocityflag = 1; + if (vzvarstr && vzvarstyle == ATOM) velocityflag = 1; + } + + maxatom = atom->nmax; + memory->destroy(displace); + memory->destroy(velocity); + if (displaceflag) memory->create(displace,maxatom,3,"move:displace"); + else displace = NULL; + if (velocityflag) memory->create(velocity,maxatom,3,"move:velocity"); + else velocity = NULL; +} + +void FixMesoMove::setup_pre_force (int /*vflag*/) { + // set vest equal to v + double **v = atom->v; + double **vest = atom->vest; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) + nlocal = atom->nfirst; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + vest[i][0] = v[i][0]; + vest[i][1] = v[i][1]; + vest[i][2] = v[i][2]; + } + } +} + +/* ---------------------------------------------------------------------- + set x,v of particles +------------------------------------------------------------------------- */ + +void FixMesoMove::initial_integrate (int /*vflag*/) { + double ddotr,dx,dy,dz; + double dtfm; + double xold[3],a[3],b[3],c[3],d[3],disp[3],disp_next[3]; + + double delta = (update->ntimestep - time_origin) * dt; + double delta_next = (update->ntimestep - time_origin + 1) * dt; + + double **x = atom->x; + double **v = atom->v; + double **vest = atom->vest; + double *rho = atom->rho; + double *drho = atom->drho; + double *e = atom->e; + double *de = atom->de; + double **f = atom->f; + double **omega = atom->omega; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int rmass_flag = atom->rmass_flag; + + if (igroup == atom->firstgroup) + nlocal = atom->nfirst; + + // for linear: X = X0 + V*dt + + if (mstyle == LINEAR) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + xold[0] = x[i][0]; + xold[1] = x[i][1]; + xold[2] = x[i][2]; + + e[i] += dtf * de[i]; // half-step update of particle internal energy + rho[i] += dtf * drho[i]; // ... and density + + if (vxflag) { + vest[i][0] = v[i][0] = vx; + x[i][0] = xoriginal[i][0] + vx*delta; + } else { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0]; + v[i][0] += dtfm * f[i][0]; + x[i][0] += dtv * v[i][0]; + } + + if (vyflag) { + vest[i][1] = v[i][1] = vy; + x[i][1] = xoriginal[i][1] + vy*delta; + } else { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1]; + v[i][1] += dtfm * f[i][1]; + x[i][1] += dtv * v[i][1]; + } + + if (vzflag) { + vest[i][2] = v[i][2] = vz; + x[i][2] = xoriginal[i][2] + vz*delta; + } else { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2]; + v[i][2] += dtfm * f[i][2]; + x[i][2] += dtv * v[i][2]; + } + + domain->remap_near(x[i],xold); + } + } + + // for wiggle: X = X0 + A sin(w*dt) + + } else if (mstyle == WIGGLE) { + double arg = omega_rotate * delta; + double arg_next = omega_rotate * delta_next; + double sine = sin(arg); + double cosine = cos(arg); + double cosine_next = cos(arg_next); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + xold[0] = x[i][0]; + xold[1] = x[i][1]; + xold[2] = x[i][2]; + + e[i] += dtf * de[i]; // half-step update of particle internal energy + rho[i] += dtf * drho[i]; // ... and density + + if (axflag) { + v[i][0] = ax*omega_rotate*cosine; + vest[i][0] = ax*omega_rotate*cosine_next; + x[i][0] = xoriginal[i][0] + ax*sine; + } else { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0]; + v[i][0] += dtfm * f[i][0]; + x[i][0] += dtv * v[i][0]; + } + + if (ayflag) { + v[i][1] = ay*omega_rotate*cosine; + vest[i][1] = ay*omega_rotate*cosine_next; + x[i][1] = xoriginal[i][1] + ay*sine; + } else { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1]; + v[i][1] += dtfm * f[i][1]; + x[i][1] += dtv * v[i][1]; + } + + if (azflag) { + v[i][2] = az*omega_rotate*cosine; + vest[i][2] = az*omega_rotate*cosine_next; + x[i][2] = xoriginal[i][2] + az*sine; + } else { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2]; + v[i][2] += dtfm * f[i][2]; + x[i][2] += dtv * v[i][2]; + } + + domain->remap_near(x[i],xold); + } + } + + // for rotate by right-hand rule around omega: + // P = point = vector = point of rotation + // R = vector = axis of rotation + // w = omega of rotation (from period) + // X0 = xoriginal = initial coord of atom + // R0 = runit = unit vector for R + // D = X0 - P = vector from P to X0 + // C = (D dot R0) R0 = projection of atom coord onto R line + // A = D - C = vector from R line to X0 + // B = R0 cross A = vector perp to A in plane of rotation + // A,B define plane of circular rotation around R line + // X = P + C + A cos(w*dt) + B sin(w*dt) + // V = w R0 cross (A cos(w*dt) + B sin(w*dt)) + + } else if (mstyle == ROTATE) { + double arg = omega_rotate * delta; + double arg_next = omega_rotate * delta_next; + double sine = sin(arg); + double cosine = cos(arg); + double sine_next = sin(arg_next); + double cosine_next = cos(arg_next); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + xold[0] = x[i][0]; + xold[1] = x[i][1]; + xold[2] = x[i][2]; + + e[i] += dtf * de[i]; // half-step update of particle internal energy + rho[i] += dtf * drho[i]; // ... and density + + d[0] = xoriginal[i][0] - point[0]; + d[1] = xoriginal[i][1] - point[1]; + d[2] = xoriginal[i][2] - point[2]; + ddotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2]; + c[0] = ddotr*runit[0]; + c[1] = ddotr*runit[1]; + c[2] = ddotr*runit[2]; + a[0] = d[0] - c[0]; + a[1] = d[1] - c[1]; + a[2] = d[2] - c[2]; + b[0] = runit[1]*a[2] - runit[2]*a[1]; + b[1] = runit[2]*a[0] - runit[0]*a[2]; + b[2] = runit[0]*a[1] - runit[1]*a[0]; + disp[0] = a[0]*cosine + b[0]*sine; + disp[1] = a[1]*cosine + b[1]*sine; + disp[2] = a[2]*cosine + b[2]*sine; + disp_next[0] = a[0]*cosine_next + b[0]*sine_next; + disp_next[1] = a[1]*cosine_next + b[1]*sine_next; + disp_next[2] = a[2]*cosine_next + b[2]*sine_next; + + x[i][0] = point[0] + c[0] + disp[0]; + x[i][1] = point[1] + c[1] + disp[1]; + x[i][2] = point[2] + c[2] + disp[2]; + v[i][0] = omega_rotate * (runit[1]*disp[2] - runit[2]*disp[1]); + v[i][1] = omega_rotate * (runit[2]*disp[0] - runit[0]*disp[2]); + v[i][2] = omega_rotate * (runit[0]*disp[1] - runit[1]*disp[0]); + vest[i][0] = omega_rotate * (runit[1]*disp_next[2] - runit[2]*disp_next[1]); + vest[i][1] = omega_rotate * (runit[2]*disp_next[0] - runit[0]*disp_next[2]); + vest[i][2] = omega_rotate * (runit[0]*disp_next[1] - runit[1]*disp_next[0]); + + domain->remap_near(x[i],xold); + } + } + + // for variable: compute x,v from variables + + } else if (mstyle == VARIABLE) { + + // reallocate displace and velocity arrays as necessary + + if ((displaceflag || velocityflag) && atom->nmax > maxatom) { + maxatom = atom->nmax; + if (displaceflag) { + memory->destroy(displace); + memory->create(displace,maxatom,3,"move:displace"); + } + if (velocityflag) { + memory->destroy(velocity); + memory->create(velocity,maxatom,3,"move:velocity"); + } + } + + // pre-compute variable values, wrap with clear/add + + modify->clearstep_compute(); + + if (xvarstr) { + if (xvarstyle == EQUAL) dx = input->variable->compute_equal(xvar); + else input->variable->compute_atom(xvar,igroup,&displace[0][0],3,0); + } + if (yvarstr) { + if (yvarstyle == EQUAL) dy = input->variable->compute_equal(yvar); + else input->variable->compute_atom(yvar,igroup,&displace[0][1],3,0); + } + if (zvarstr) { + if (zvarstyle == EQUAL) dz = input->variable->compute_equal(zvar); + else input->variable->compute_atom(zvar,igroup,&displace[0][2],3,0); + } + if (vxvarstr) { + if (vxvarstyle == EQUAL) vx = input->variable->compute_equal(vxvar); + else input->variable->compute_atom(vxvar,igroup,&velocity[0][0],3,0); + } + if (vyvarstr) { + if (vyvarstyle == EQUAL) vy = input->variable->compute_equal(vyvar); + else input->variable->compute_atom(vyvar,igroup,&velocity[0][1],3,0); + } + if (vzvarstr) { + if (vzvarstyle == EQUAL) vz = input->variable->compute_equal(vzvar); + else input->variable->compute_atom(vzvar,igroup,&velocity[0][2],3,0); + } + + modify->addstep_compute(update->ntimestep + 1); + + // update x,v + // vest (velocity in next step) could be different from v in the next + // step, but this is the best we could do + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + xold[0] = x[i][0]; + xold[1] = x[i][1]; + xold[2] = x[i][2]; + + if (xvarstr && vxvarstr) { + if (vxvarstyle == EQUAL) { + vest[i][0] = 2*vx - v[i][0]; + v[i][0] = vx; + } + else { + vest[i][0] = 2*velocity[i][0] - v[i][0]; + v[i][0] = velocity[i][0]; + } + if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx; + else x[i][0] = xoriginal[i][0] + displace[i][0]; + } else if (xvarstr) { + if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx; + else x[i][0] = xoriginal[i][0] + displace[i][0]; + } else if (vxvarstr) { + if (vxvarstyle == EQUAL) { + vest[i][0] = 2*vx - v[i][0]; + v[i][0] = vx; + } + else { + vest[i][0] = 2*velocity[i][0] - v[i][0]; + v[i][0] = velocity[i][0]; + } + x[i][0] += dtv * v[i][0]; + } else { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0]; + v[i][0] += dtfm * f[i][0]; + x[i][0] += dtv * v[i][0]; + } + + if (yvarstr && vyvarstr) { + if (vyvarstyle == EQUAL) { + vest[i][1] = 2*vy - v[i][1]; + v[i][1] = vy; + } + else { + vest[i][1] = 2*velocity[i][1] - v[i][1]; + v[i][1] = velocity[i][1]; + } + if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy; + else x[i][1] = xoriginal[i][1] + displace[i][1]; + } else if (yvarstr) { + if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy; + else x[i][1] = xoriginal[i][1] + displace[i][1]; + } else if (vyvarstr) { + if (vyvarstyle == EQUAL) { + vest[i][1] = 2*vy - v[i][1]; + v[i][1] = vy; + } + else { + vest[i][1] = 2*velocity[i][1] - v[i][1]; + v[i][1] = velocity[i][1]; + } + x[i][1] += dtv * v[i][1]; + } else { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1]; + v[i][1] += dtfm * f[i][1]; + x[i][1] += dtv * v[i][1]; + } + + if (zvarstr && vzvarstr) { + if (vzvarstyle == EQUAL) { + vest[i][2] = 2*vz - v[i][2]; + v[i][2] = vz; + } + else { + vest[i][2] = 2*velocity[i][2] - v[i][2]; + v[i][2] = velocity[i][2]; + } + if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz; + else x[i][2] = xoriginal[i][2] + displace[i][2]; + } else if (zvarstr) { + if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz; + else x[i][2] = xoriginal[i][2] + displace[i][2]; + } else if (vzvarstr) { + if (vzvarstyle == EQUAL) { + vest[i][2] = 2*vz - v[i][2]; + v[i][2] = vz; + } + else { + vest[i][2] = 2*velocity[i][2] - v[i][2]; + v[i][2] = velocity[i][2]; + } + x[i][2] += dtv * v[i][2]; + } else { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2]; + v[i][2] += dtfm * f[i][2]; + x[i][2] += dtv * v[i][2]; + } + + domain->remap_near(x[i],xold); + } + } + } +} + +/* ---------------------------------------------------------------------- + final NVE of particles with NULL components +------------------------------------------------------------------------- */ + +void FixMesoMove::final_integrate () { + double dtfm; + + int xflag = 1; + if (mstyle == LINEAR && vxflag) xflag = 0; + else if (mstyle == WIGGLE && axflag) xflag = 0; + else if (mstyle == ROTATE) xflag = 0; + else if (mstyle == VARIABLE && (xvarstr || vxvarstr)) xflag = 0; + + int yflag = 1; + if (mstyle == LINEAR && vyflag) yflag = 0; + else if (mstyle == WIGGLE && ayflag) yflag = 0; + else if (mstyle == ROTATE) yflag = 0; + else if (mstyle == VARIABLE && (yvarstr || vyvarstr)) yflag = 0; + + int zflag = 1; + if (mstyle == LINEAR && vzflag) zflag = 0; + else if (mstyle == WIGGLE && azflag) zflag = 0; + else if (mstyle == ROTATE) zflag = 0; + else if (mstyle == VARIABLE && (zvarstr || vzvarstr)) zflag = 0; + + double **v = atom->v; + double **f = atom->f; + double *e = atom->e; + double *de = atom->de; + double *rho = atom->rho; + double *drho = atom->drho; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int rmass_flag = atom->rmass_flag; + + if (igroup == atom->firstgroup) + nlocal = atom->nfirst; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + e[i] += dtf * de[i]; + rho[i] += dtf * drho[i]; + + if (xflag) { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + v[i][0] += dtfm * f[i][0]; + } + + if (yflag) { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + v[i][1] += dtfm * f[i][1]; + } + + if (zflag) { + dtfm = rmass_flag ? dtf / rmass[i] : dtf / mass[type[i]]; + v[i][2] += dtfm * f[i][2]; + } + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double FixMesoMove::memory_usage () { + double bytes = atom->nmax*3 * sizeof(double); + if (displaceflag) bytes += atom->nmax*3 * sizeof(double); + if (velocityflag) bytes += atom->nmax*3 * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + pack entire state of Fix into one write +------------------------------------------------------------------------- */ + +void FixMesoMove::write_restart (FILE *fp) { + int n = 0; + double list[1]; + list[n++] = time_origin; + + if (comm->me == 0) { + int size = n * sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(list,sizeof(double),n,fp); + } +} + +/* ---------------------------------------------------------------------- + use state info from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixMesoMove::restart (char *buf) { + int n = 0; + double *list = (double *) buf; + + time_origin = static_cast (list[n++]); +} + +/* ---------------------------------------------------------------------- + allocate atom-based array +------------------------------------------------------------------------- */ + +void FixMesoMove::grow_arrays (int nmax) { + memory->grow(xoriginal,nmax,3,"move:xoriginal"); + array_atom = xoriginal; +} + +/* ---------------------------------------------------------------------- + copy values within local atom-based array +------------------------------------------------------------------------- */ + +void FixMesoMove::copy_arrays (int i, int j, int /*delflag*/) { + xoriginal[j][0] = xoriginal[i][0]; + xoriginal[j][1] = xoriginal[i][1]; + xoriginal[j][2] = xoriginal[i][2]; +} + +/* ---------------------------------------------------------------------- + initialize one atom's array values, called when atom is created +------------------------------------------------------------------------- */ + +void FixMesoMove::set_arrays (int i) { + double **x = atom->x; + imageint *image = atom->image; + int *mask = atom->mask; + + // particle not in group + + if (!(mask[i] & groupbit)) { + xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; + return; + } + + // current time still equal fix creation time + + if (update->ntimestep == time_origin) { + domain->unmap(x[i],image[i],xoriginal[i]); + return; + } + + // backup particle to time_origin + + if (mstyle == VARIABLE) + error->all(FLERR,"Cannot add atoms to fix meso/move variable"); + + domain->unmap(x[i],image[i],xoriginal[i]); + double delta = (update->ntimestep - time_origin) * update->dt; + + if (mstyle == LINEAR) { + if (vxflag) xoriginal[i][0] -= vx * delta; + if (vyflag) xoriginal[i][1] -= vy * delta; + if (vzflag) xoriginal[i][2] -= vz * delta; + } else if (mstyle == WIGGLE) { + double arg = omega_rotate * delta; + double sine = sin(arg); + if (axflag) xoriginal[i][0] -= ax*sine; + if (ayflag) xoriginal[i][1] -= ay*sine; + if (azflag) xoriginal[i][2] -= az*sine; + } else if (mstyle == ROTATE) { + double a[3],b[3],c[3],d[3],disp[3],ddotr; + double arg = - omega_rotate * delta; + double sine = sin(arg); + double cosine = cos(arg); + d[0] = x[i][0] - point[0]; + d[1] = x[i][1] - point[1]; + d[2] = x[i][2] - point[2]; + ddotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2]; + c[0] = ddotr*runit[0]; + c[1] = ddotr*runit[1]; + c[2] = ddotr*runit[2]; + + a[0] = d[0] - c[0]; + a[1] = d[1] - c[1]; + a[2] = d[2] - c[2]; + b[0] = runit[1]*a[2] - runit[2]*a[1]; + b[1] = runit[2]*a[0] - runit[0]*a[2]; + b[2] = runit[0]*a[1] - runit[1]*a[0]; + disp[0] = a[0]*cosine + b[0]*sine; + disp[1] = a[1]*cosine + b[1]*sine; + disp[2] = a[2]*cosine + b[2]*sine; + + xoriginal[i][0] = point[0] + c[0] + disp[0]; + xoriginal[i][1] = point[1] + c[1] + disp[1]; + xoriginal[i][2] = point[2] + c[2] + disp[2]; + } +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based array for exchange with another proc +------------------------------------------------------------------------- */ + +int FixMesoMove::pack_exchange (int i, double *buf) { + buf[0] = xoriginal[i][0]; + buf[1] = xoriginal[i][1]; + buf[2] = xoriginal[i][2]; + return 3; +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based array from exchange with another proc +------------------------------------------------------------------------- */ + +int FixMesoMove::unpack_exchange (int nlocal, double *buf) { + xoriginal[nlocal][0] = buf[0]; + xoriginal[nlocal][1] = buf[1]; + xoriginal[nlocal][2] = buf[2]; + return 3; +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based arrays for restart file +------------------------------------------------------------------------- */ + +int FixMesoMove::pack_restart (int i, double *buf) { + buf[0] = 4; + buf[1] = xoriginal[i][0]; + buf[2] = xoriginal[i][1]; + buf[3] = xoriginal[i][2]; + return 4; +} + +/* ---------------------------------------------------------------------- + unpack values from atom->extra array to restart the fix +------------------------------------------------------------------------- */ + +void FixMesoMove::unpack_restart (int nlocal, int nth) { + double **extra = atom->extra; + + // skip to Nth set of extra values + + int m = 0; + for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); + m++; + + xoriginal[nlocal][0] = extra[nlocal][m++]; + xoriginal[nlocal][1] = extra[nlocal][m++]; + xoriginal[nlocal][2] = extra[nlocal][m++]; +} + +/* ---------------------------------------------------------------------- + maxsize of any atom's restart data +------------------------------------------------------------------------- */ + +int FixMesoMove::maxsize_restart () { + return 4; +} + +/* ---------------------------------------------------------------------- + size of atom nlocal's restart data +------------------------------------------------------------------------- */ + +int FixMesoMove::size_restart (int nlocal) { + return 4; +} + +/* ---------------------------------------------------------------------- */ + +void FixMesoMove::reset_dt () { + error->all(FLERR,"Resetting timestep size is not allowed with fix meso/move"); +} diff --git a/src/USER-SDPD/fix_meso_move.h b/src/USER-SDPD/fix_meso_move.h new file mode 100644 index 0000000000..a5e213ba81 --- /dev/null +++ b/src/USER-SDPD/fix_meso_move.h @@ -0,0 +1,127 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(meso/move,FixMesoMove) + +#else + +#ifndef LMP_FIX_MESO_MOVE_H +#define LMP_FIX_MESO_MOVE_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixMesoMove : public Fix { + public: + FixMesoMove (class LAMMPS *, int, char **); + ~FixMesoMove (); + int setmask (); + void init (); + void setup_pre_force (int); + void initial_integrate (int); + void final_integrate (); + + double memory_usage (); + void write_restart (FILE *); + void restart (char *); + void grow_arrays (int); + void copy_arrays (int, int, int); + void set_arrays (int); + int pack_exchange (int, double *); + int unpack_exchange (int, double *); + int pack_restart (int, double *); + void unpack_restart (int, int); + int maxsize_restart (); + int size_restart (int); + + void reset_dt (); + + private: + char *xvarstr,*yvarstr,*zvarstr,*vxvarstr,*vyvarstr,*vzvarstr; + int mstyle; + int vxflag,vyflag,vzflag,axflag,ayflag,azflag; + double vx,vy,vz,ax,ay,az; + double period,omega_rotate; + double point[3],axis[3],runit[3]; + double dt,dtv,dtf; + int xvar,yvar,zvar,vxvar,vyvar,vzvar; + int xvarstyle,yvarstyle,zvarstyle,vxvarstyle,vyvarstyle,vzvarstyle; + int time_origin; + + double **xoriginal; // original coords of atoms + int displaceflag,velocityflag; + int maxatom; + double **displace,**velocity; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Fix meso/move cannot set linear z motion for 2d problem + +Self-explanatory. + +E: Fix meso/move cannot set wiggle z motion for 2d problem + +Self-explanatory. + +E: Fix meso/move cannot rotate aroung non z-axis for 2d problem + +Self-explanatory. + +E: Fix meso/move cannot define z or vz variable for 2d problem + +Self-explanatory. + +W: Fix meso/move does not update angular momentum + +Atoms store this quantity, but fix meso/move does not (yet) update it. + +W: Fix meso/move does not update quaternions + +Atoms store this quantity, but fix meso/move does not (yet) update it. + +E: Zero length rotation vector with fix meso/move + +Self-explanatory. + +E: Variable name for fix meso/move does not exist + +Self-explanatory. + +E: Variable for fix meso/move is invalid style + +Only equal-style variables can be used. + +E: Cannot add atoms to fix meso/move variable + +Atoms can not be added afterwards to this fix option. + +E: Resetting timestep size is not allowed with fix meso/move + +This is because fix meso/move is moving atoms based on elapsed time. + +*/ diff --git a/src/USER-SDPD/fix_rigid_meso.cpp b/src/USER-SDPD/fix_rigid_meso.cpp new file mode 100644 index 0000000000..fd881852f7 --- /dev/null +++ b/src/USER-SDPD/fix_rigid_meso.cpp @@ -0,0 +1,500 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan) + references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) + Miller et al., J Chem Phys. 116, 8649-8659 (2002) +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: + Morteza Jalalvand (IASBS) jalalvand.m AT gmail.com + + This is an extension of fix/rigid/nve to SPH/SDPD particles + You can see the original copyright notice of fix/rigid authors above + Note that the Kamberaj paper was related to the nvt variant + and all codes relevant to that has been removed +------------------------------------------------------------------------- */ + +#include +#include "fix_rigid_meso.h" +#include "math_extra.h" +#include "atom.h" +#include "compute.h" +#include "domain.h" +#include "update.h" +#include "modify.h" +#include "group.h" +#include "force.h" +#include "output.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixRigidMeso::FixRigidMeso (LAMMPS *lmp, int narg, char **arg) : +FixRigid (lmp, narg, arg) { + scalar_flag = 0; + size_array_cols = 28; + if ((atom->e_flag != 1) || (atom->rho_flag != 1)) + error->all (FLERR, "fix rigid/meso command requires atom_style with" + " both energy and density"); + + if (langflag || tstat_flag) + error->all (FLERR,"Can not use thermostat with fix rigid/meso"); + + if (pstat_flag) + error->all (FLERR,"Can not use barostat with fix rigid/meso"); + + // memory allocation and initialization + + memory->create(conjqm,nbody,4,"rigid_nh:conjqm"); +} + +/* ---------------------------------------------------------------------- */ + +FixRigidMeso::~FixRigidMeso () { + memory->destroy(conjqm); +} + +/* ---------------------------------------------------------------------- */ + +int FixRigidMeso::setmask () { + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + mask |= PRE_NEIGHBOR; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidMeso::setup (int vflag) { + FixRigid::setup(vflag); + + double mbody[3]; + for (int ibody = 0; ibody < nbody; ibody++) { + MathExtra::transpose_matvec (ex_space[ibody],ey_space[ibody],ez_space[ibody], + angmom[ibody],mbody); + MathExtra::quatvec (quat[ibody],mbody,conjqm[ibody]); + conjqm[ibody][0] *= 2.0; + conjqm[ibody][1] *= 2.0; + conjqm[ibody][2] *= 2.0; + conjqm[ibody][3] *= 2.0; + } +} + +/* ---------------------------------------------------------------------- + perform preforce velocity Verlet integration + see Kamberaj paper for step references +------------------------------------------------------------------------- */ + +void FixRigidMeso::initial_integrate (int vflag) { + double dtfm,mbody[3],tbody[3],fquat[4]; + double dtf2 = dtf * 2.0; + + // update xcm, vcm, quat, conjqm and angmom + + for (int ibody = 0; ibody < nbody; ibody++) { + + // step 1.1 - update vcm by 1/2 step + + dtfm = dtf / masstotal[ibody]; + vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; + vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; + vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; + + // step 1.2 - update xcm by full step + + xcm[ibody][0] += dtv * vcm[ibody][0]; + xcm[ibody][1] += dtv * vcm[ibody][1]; + xcm[ibody][2] += dtv * vcm[ibody][2]; + + // step 1.3 - apply torque (body coords) to quaternion momentum + + torque[ibody][0] *= tflag[ibody][0]; + torque[ibody][1] *= tflag[ibody][1]; + torque[ibody][2] *= tflag[ibody][2]; + + MathExtra::transpose_matvec (ex_space[ibody],ey_space[ibody],ez_space[ibody], + torque[ibody],tbody); + MathExtra::quatvec (quat[ibody],tbody,fquat); + + conjqm[ibody][0] += dtf2 * fquat[0]; + conjqm[ibody][1] += dtf2 * fquat[1]; + conjqm[ibody][2] += dtf2 * fquat[2]; + conjqm[ibody][3] += dtf2 * fquat[3]; + + // step 1.4 to 1.13 - use no_squish rotate to update p and q + + MathExtra::no_squish_rotate (3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + MathExtra::no_squish_rotate (2,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + MathExtra::no_squish_rotate (1,conjqm[ibody],quat[ibody],inertia[ibody],dtv); + MathExtra::no_squish_rotate (2,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + MathExtra::no_squish_rotate (3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + + // update exyz_space + // transform p back to angmom + // update angular velocity + + MathExtra::q_to_exyz (quat[ibody],ex_space[ibody],ey_space[ibody], + ez_space[ibody]); + MathExtra::invquatvec (quat[ibody],conjqm[ibody],mbody); + MathExtra::matvec (ex_space[ibody],ey_space[ibody],ez_space[ibody], + mbody,angmom[ibody]); + + angmom[ibody][0] *= 0.5; + angmom[ibody][1] *= 0.5; + angmom[ibody][2] *= 0.5; + + MathExtra::angmom_to_omega (angmom[ibody],ex_space[ibody],ey_space[ibody], + ez_space[ibody],inertia[ibody],omega[ibody]); + } + + // virial setup before call to set_xv + + if (vflag) v_setup(vflag); + else evflag = 0; + + // set coords/orient and velocity/rotation of atoms in rigid bodies + // from quarternion and omega + + set_xv(); +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidMeso::final_integrate () { + int i,ibody; + double dtfm,xy,xz,yz; + double mbody[3],tbody[3],fquat[4]; + + double dtf2 = dtf * 2.0; + + // late calculation of forces and torques (if requested) + + if (!earlyflag) compute_forces_and_torques(); + + // update vcm and angmom + // fflag,tflag = 0 for some dimensions in 2d + + for (ibody = 0; ibody < nbody; ibody++) { + + // update vcm by 1/2 step + + dtfm = dtf / masstotal[ibody]; + + vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; + vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; + vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; + + // update conjqm, then transform to angmom, set velocity again + // virial is already setup from initial_integrate + + torque[ibody][0] *= tflag[ibody][0]; + torque[ibody][1] *= tflag[ibody][1]; + torque[ibody][2] *= tflag[ibody][2]; + + MathExtra::transpose_matvec (ex_space[ibody],ey_space[ibody], + ez_space[ibody],torque[ibody],tbody); + MathExtra::quatvec (quat[ibody],tbody,fquat); + + conjqm[ibody][0] += dtf2 * fquat[0]; + conjqm[ibody][1] += dtf2 * fquat[1]; + conjqm[ibody][2] += dtf2 * fquat[2]; + conjqm[ibody][3] += dtf2 * fquat[3]; + + MathExtra::invquatvec (quat[ibody],conjqm[ibody],mbody); + MathExtra::matvec (ex_space[ibody],ey_space[ibody],ez_space[ibody], + mbody,angmom[ibody]); + + angmom[ibody][0] *= 0.5; + angmom[ibody][1] *= 0.5; + angmom[ibody][2] *= 0.5; + + MathExtra::angmom_to_omega (angmom[ibody],ex_space[ibody],ey_space[ibody], + ez_space[ibody],inertia[ibody],omega[ibody]); + } + + // set velocity/rotation of atoms in rigid bodies + // virial is already setup from initial_integrate + + set_v(); +} + +/* ---------------------------------------------------------------------- + set space-frame coords and velocity of each atom in each rigid body + set orientation and rotation of extended particles + x = Q displace + Xcm, mapped back to periodic box + v = Vcm + (W cross (x - Xcm)) +------------------------------------------------------------------------- */ + +void FixRigidMeso::set_xv () { + int ibody; + int xbox,ybox,zbox; + double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; + double xy,xz,yz; + double ione[3],exone[3],eyone[3],ezone[3],vr[6],p[3][3]; + + double **x = atom->x; + double **v = atom->v; + double **vest = atom->vest; + double **f = atom->f; + double *e = atom->e; + double *de = atom->de; + double *rho = atom->rho; + double *drho = atom->drho; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int nlocal = atom->nlocal; + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + + if (triclinic) { + xy = domain->xy; + xz = domain->xz; + yz = domain->yz; + } + + // set x and v of each atom + + for (int i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + + // half-step update of particle internal energy and density + e[i] += dtf * de[i]; + rho[i] += dtf * drho[i]; + + ibody = body[i]; + + xbox = (xcmimage[i] & IMGMASK) - IMGMAX; + ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX; + + // save old positions and velocities for virial + + if (evflag) { + if (triclinic == 0) { + x0 = x[i][0] + xbox*xprd; + x1 = x[i][1] + ybox*yprd; + x2 = x[i][2] + zbox*zprd; + } else { + x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; + x1 = x[i][1] + ybox*yprd + zbox*yz; + x2 = x[i][2] + zbox*zprd; + } + } + + v0 = v[i][0]; + v1 = v[i][1]; + v2 = v[i][2]; + + // x = displacement from center-of-mass, based on body orientation + // v = vcm + omega around center-of-mass + // vest = 2*v - v_old + + MathExtra::matvec (ex_space[ibody],ey_space[ibody], + ez_space[ibody],displace[i],x[i]); + + v[i][0] = omega[ibody][1]*x[i][2] - omega[ibody][2]*x[i][1] + + vcm[ibody][0]; + v[i][1] = omega[ibody][2]*x[i][0] - omega[ibody][0]*x[i][2] + + vcm[ibody][1]; + v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] + + vcm[ibody][2]; + + vest[i][0] = 2*v[i][0] - v0; + vest[i][1] = 2*v[i][1] - v1; + vest[i][2] = 2*v[i][2] - v2; + + // add center of mass to displacement + // map back into periodic box via xbox,ybox,zbox + // for triclinic, add in box tilt factors as well + + if (triclinic == 0) { + x[i][0] += xcm[ibody][0] - xbox*xprd; + x[i][1] += xcm[ibody][1] - ybox*yprd; + x[i][2] += xcm[ibody][2] - zbox*zprd; + } else { + x[i][0] += xcm[ibody][0] - xbox*xprd - ybox*xy - zbox*xz; + x[i][1] += xcm[ibody][1] - ybox*yprd - zbox*yz; + x[i][2] += xcm[ibody][2] - zbox*zprd; + } + + // virial = unwrapped coords dotted into body constraint force + // body constraint force = implied force due to v change minus f external + // assume f does not include forces internal to body + // 1/2 factor b/c final_integrate contributes other half + // assume per-atom contribution is due to constraint force on that atom + + if (evflag) { + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; + fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; + fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; + + vr[0] = 0.5*x0*fc0; + vr[1] = 0.5*x1*fc1; + vr[2] = 0.5*x2*fc2; + vr[3] = 0.5*x0*fc1; + vr[4] = 0.5*x0*fc2; + vr[5] = 0.5*x1*fc2; + + v_tally(1,&i,1.0,vr); + } + } + + // set orientation, omega, angmom of each extended particle + + if (extended) { + // TBD + } +} + +/* ---------------------------------------------------------------------- + set space-frame velocity of each atom in a rigid body + set omega and angmom of extended particles + v = Vcm + (W cross (x - Xcm)) +------------------------------------------------------------------------- */ + +void FixRigidMeso::set_v () { + int xbox,ybox,zbox; + double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; + double xy,xz,yz; + double ione[3],exone[3],eyone[3],ezone[3],delta[3],vr[6]; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *e = atom->e; + double *de = atom->de; + double *rho = atom->rho; + double *drho = atom->drho; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int nlocal = atom->nlocal; + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + if (triclinic) { + xy = domain->xy; + xz = domain->xz; + yz = domain->yz; + } + + // set v of each atom + + for (int i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + + // half-step update of particle internal energy and density + e[i] += dtf * de[i]; + rho[i] += dtf * drho[i]; + + const int ibody = body[i]; + + MathExtra::matvec (ex_space[ibody],ey_space[ibody], + ez_space[ibody],displace[i],delta); + + // save old velocities for virial + + if (evflag) { + v0 = v[i][0]; + v1 = v[i][1]; + v2 = v[i][2]; + } + + v[i][0] = omega[ibody][1]*delta[2] - omega[ibody][2]*delta[1] + + vcm[ibody][0]; + v[i][1] = omega[ibody][2]*delta[0] - omega[ibody][0]*delta[2] + + vcm[ibody][1]; + v[i][2] = omega[ibody][0]*delta[1] - omega[ibody][1]*delta[0] + + vcm[ibody][2]; + + // virial = unwrapped coords dotted into body constraint force + // body constraint force = implied force due to v change minus f external + // assume f does not include forces internal to body + // 1/2 factor b/c initial_integrate contributes other half + // assume per-atom contribution is due to constraint force on that atom + + if (evflag) { + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; + fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; + fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; + + xbox = (xcmimage[i] & IMGMASK) - IMGMAX; + ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX; + + if (triclinic == 0) { + x0 = x[i][0] + xbox*xprd; + x1 = x[i][1] + ybox*yprd; + x2 = x[i][2] + zbox*zprd; + } else { + x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; + x1 = x[i][1] + ybox*yprd + zbox*yz; + x2 = x[i][2] + zbox*zprd; + } + + vr[0] = 0.5*x0*fc0; + vr[1] = 0.5*x1*fc1; + vr[2] = 0.5*x2*fc2; + vr[3] = 0.5*x0*fc1; + vr[4] = 0.5*x0*fc2; + vr[5] = 0.5*x1*fc2; + + v_tally(1,&i,1.0,vr); + } + } + + // set omega, angmom of each extended particle + + if (extended) { + // TBD + } +} + +/* ---------------------------------------------------------------------- + return attributes of a rigid body + 19 values per body + xcm = 0,1,2; vcm = 3,4,5; fcm = 6,7,8; + quat = 9,10,11,12; omega = 13,14,15; torque = 16,17,18; + inertia = 19,20,21; angmom = 22,23,24; + image = 25,26,27 +------------------------------------------------------------------------- */ + +double FixRigidMeso::compute_array (int i, int j) { + if (j < 3) return xcm[i][j]; + if (j < 6) return vcm[i][j-3]; + if (j < 9) return fcm[i][j-6]; + if (j < 13) return quat[i][j-9]; + if (j < 16) return omega[i][j-13]; + if (j < 19) return torque[i][j-16]; + if (j < 22) return inertia[i][j-19]; + if (j < 25) return angmom[i][j-22]; + if (j == 25) return (imagebody[i] & IMGMASK) - IMGMAX; + if (j == 26) return (imagebody[i] >> IMGBITS & IMGMASK) - IMGMAX; + return (imagebody[i] >> IMG2BITS) - IMGMAX; +} diff --git a/src/USER-SDPD/fix_rigid_meso.h b/src/USER-SDPD/fix_rigid_meso.h new file mode 100644 index 0000000000..ac73a9503b --- /dev/null +++ b/src/USER-SDPD/fix_rigid_meso.h @@ -0,0 +1,69 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(rigid/meso,FixRigidMeso) + +#else + +#ifndef LMP_FIX_RIGID_MESO_H +#define LMP_FIX_RIGID_MESO_H + +#include "fix_rigid.h" + +namespace LAMMPS_NS { + +class FixRigidMeso : public FixRigid { + public: + FixRigidMeso (class LAMMPS *, int, char **); + ~FixRigidMeso (); + int setmask (); + void setup (int); + void initial_integrate (int); + void final_integrate (); + double compute_scalar () {} + double compute_array (int, int); + + protected: + void set_xv (); + void set_v (); + double **conjqm; // conjugate quaternion momentum +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: fix rigid/meso command requires atom_style with both energy and density + +You should use atom_style meso with this fix + +E: Can not use thermostat with fix rigid/meso + +Self-explanatory + +E: Can not use barostat with fix rigid/meso + +Self-explanatory + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +*/ diff --git a/src/USER-SDPD/pair_sdpd_taitwater_isothermal.cpp b/src/USER-SDPD/pair_sdpd_taitwater_isothermal.cpp new file mode 100644 index 0000000000..d8c32e7f6c --- /dev/null +++ b/src/USER-SDPD/pair_sdpd_taitwater_isothermal.cpp @@ -0,0 +1,321 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: + Morteza Jalalvand (IASBS) jalalvand.m AT gmail.com + + references: Espanol and Revenga, Phys Rev E 67, 026705 (2003) +------------------------------------------------------------------------- */ + +#include +#include +#include "pair_sdpd_taitwater_isothermal.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" +#include "domain.h" +#include "update.h" +#ifndef USE_ZEST +#include "random_mars.h" +#endif + +using namespace LAMMPS_NS; + +static const double sqrt_2_inv = std::sqrt(0.5); + +/* ---------------------------------------------------------------------- */ + +PairSDPDTaitwaterIsothermal::PairSDPDTaitwaterIsothermal (LAMMPS *lmp) +: Pair (lmp) { + restartinfo = 0; +} + +/* ---------------------------------------------------------------------- */ + +PairSDPDTaitwaterIsothermal::~PairSDPDTaitwaterIsothermal () { + if (allocated) { + memory->destroy (setflag); + memory->destroy (cutsq); + + memory->destroy (cut); + memory->destroy (rho0); + memory->destroy (soundspeed); + memory->destroy (B); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairSDPDTaitwaterIsothermal::compute (int eflag, int vflag) { + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + + int *ilist, *jlist, *numneigh, **firstneigh; + double vxtmp, vytmp, vztmp, imass, jmass, fi, fj, fvisc; + double h, ih, ihsq, velx, vely, velz; + double rsq, tmp, wfd, delVdotDelR, deltaE; + double prefactor, wiener[3][3], f_random[3]; + + if (eflag || vflag) ev_setup (eflag, vflag); + else evflag = vflag_fdotr = 0; + + double **v = atom->vest; + double **x = atom->x; + double **f = atom->f; + double *rho = atom->rho; + double *mass = atom->mass; + double *de = atom->de; + double *drho = atom->drho; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + int dimension = domain->dimension; + double dtinv = 1.0 / update->dt; + double kBoltzmann = force->boltz; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + vxtmp = v[i][0]; + vytmp = v[i][1]; + vztmp = v[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + imass = mass[itype]; + + // compute pressure of atom i with Tait EOS + tmp = rho[i] / rho0[itype]; + fi = tmp * tmp * tmp; + fi = B[itype] * (fi * fi * tmp - 1.0) / (rho[i] * rho[i]); + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + jmass = mass[jtype]; + + if (rsq < cutsq[itype][jtype]) { + h = cut[itype][jtype]; + ih = 1.0 / h; + ihsq = ih * ih; + + double r = sqrt (rsq); + wfd = h - r; + if (dimension == 3) { + // Lucy Kernel, 3d + // Note that wfd, the derivative of the weight function with respect to r, + // is lacking a factor of r. + // The missing factor of r is recovered by + // (1) using delV . delX instead of delV . (delX/r) and + // (2) using f[i][0] += delx * fpair instead of f[i][0] += (delx/r) * fpair + wfd = -25.066903536973515383e0 * wfd * wfd * ihsq * ihsq * ihsq * ih; + } else { + // Lucy Kernel, 2d + wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq; + } + + // compute pressure of atom j with Tait EOS + tmp = rho[j] / rho0[jtype]; + fj = tmp * tmp * tmp; + fj = B[jtype] * (fj * fj * tmp - 1.0) / (rho[j] * rho[j]); + + velx=vxtmp - v[j][0]; + vely=vytmp - v[j][1]; + velz=vztmp - v[j][2]; + + // dot product of velocity delta and distance vector + delVdotDelR = delx * velx + dely * vely + delz * velz; + + // Espanol Viscosity (Espanol, 2003) + + fvisc = (5. / 3.) * viscosity * imass * jmass * wfd / (rho[i]*rho[j]); + + // total pair force + fpair = -imass * jmass * (fi + fj) * wfd; + + // random force calculation + // independent increments of a Wiener process matrix +#ifdef USE_ZEST + wiener[0][0] = gaussian (generator); + wiener[1][1] = gaussian (generator); + wiener[2][2] = gaussian (generator); + + wiener[0][1] = wiener[1][0] = sqrt_2_inv * gaussian (generator); + wiener[0][2] = wiener[2][0] = sqrt_2_inv * gaussian (generator); + wiener[1][2] = wiener[2][1] = sqrt_2_inv * gaussian (generator); +#else + wiener[0][0] = random->gaussian (); + wiener[1][1] = random->gaussian (); + wiener[2][2] = random->gaussian (); + + wiener[0][1] = wiener[1][0] = sqrt_2_inv * random->gaussian (); + wiener[0][2] = wiener[2][0] = sqrt_2_inv * random->gaussian (); + wiener[1][2] = wiener[2][1] = sqrt_2_inv * random->gaussian (); +#endif + + prefactor = sqrt (-4. * kBoltzmann*temperature * fvisc * dtinv) / r; + + f_random[0] = prefactor * (wiener[0][0]*delx + wiener[0][1]*dely + wiener[0][2]*delz); + f_random[1] = prefactor * (wiener[1][0]*delx + wiener[1][1]*dely + wiener[1][2]*delz); + f_random[2] = prefactor * (wiener[2][0]*delx + wiener[2][1]*dely + wiener[2][2]*delz); + + f[i][0] += delx * fpair + (velx + delx * delVdotDelR / rsq) * fvisc + f_random[0]; + f[i][1] += dely * fpair + (vely + dely * delVdotDelR / rsq) * fvisc + f_random[1]; + f[i][2] += delz * fpair + (velz + delz * delVdotDelR / rsq) * fvisc + f_random[2]; + + // and change in density + drho[i] += jmass * delVdotDelR * wfd; + + if (newton_pair || j < nlocal) { + f[j][0] -= delx * fpair + (velx + delx * delVdotDelR / rsq) * fvisc + f_random[0]; + f[j][1] -= dely * fpair + (vely + dely * delVdotDelR / rsq) * fvisc + f_random[1]; + f[j][2] -= delz * fpair + (velz + delz * delVdotDelR / rsq) * fvisc + f_random[2]; + drho[j] += imass * delVdotDelR * wfd; + } + + if (evflag) + ev_tally (i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute (); +} + +/* ---------------------------------------------------------------------- + allocate all arrays + ------------------------------------------------------------------------- */ + +void PairSDPDTaitwaterIsothermal::allocate () { + allocated = 1; + int n = atom->ntypes; + + memory->create (setflag, n + 1, n + 1, "pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create (cutsq, n + 1, n + 1, "pair:cutsq"); + + memory->create (rho0, n + 1, "pair:rho0"); + memory->create (soundspeed, n + 1, "pair:soundspeed"); + memory->create (B, n + 1, "pair:B"); + memory->create (cut, n + 1, n + 1, "pair:cut"); +} + +/* ---------------------------------------------------------------------- + global settings + ------------------------------------------------------------------------- */ + +void PairSDPDTaitwaterIsothermal::settings (int narg, char **arg) { + if (narg != 2 && narg != 3) error->all (FLERR, "Illegal number of setting " + "arguments for pair_style sdpd/taitwater/morris/isothermal"); + + temperature = force->numeric (FLERR, arg[0]); + viscosity = force->numeric (FLERR, arg[1]); + + if (temperature <= 0) error->all (FLERR, "Temperature must be positive"); + if (viscosity <= 0) error->all (FLERR, "Viscosity must be positive"); + + // seed is immune to underflow/overflow because it is unsigned + seed = comm->nprocs + comm->me + atom->nlocal; + if (narg == 3) seed += force->inumeric (FLERR, arg[2]); +#ifdef USE_ZEST + generator.seed (seed); +#else + random = new RanMars (lmp, seed); +#endif +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs + ------------------------------------------------------------------------- */ + +void PairSDPDTaitwaterIsothermal::coeff (int narg, char **arg) { + if (narg != 5) error->all (FLERR, "Incorrect args for pair_style " + "sph/taitwater/morris coefficients"); + + if (!allocated) allocate(); + + int ilo, ihi, jlo, jhi; + force->bounds (FLERR, arg[0], atom->ntypes, ilo, ihi); + force->bounds (FLERR, arg[1], atom->ntypes, jlo, jhi); + + double rho0_one = force->numeric (FLERR,arg[2]); + double soundspeed_one = force->numeric (FLERR,arg[3]); + double cut_one = force->numeric (FLERR,arg[4]); + double B_one = soundspeed_one * soundspeed_one * rho0_one / 7.0; + + if (rho0_one <= 0) error->all (FLERR, "Density must be positive"); + if (soundspeed_one <= 0) error->all (FLERR, "Sound speed must be positive"); + if (cut_one <= 0) error->all (FLERR, "Cutoff must be positive"); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + rho0[i] = rho0_one; + soundspeed[i] = soundspeed_one; + B[i] = B_one; + for (int j = MAX(jlo,i); j <= jhi; j++) { + cut[i][j] = cut_one; + + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) + error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i + ------------------------------------------------------------------------- */ + +double PairSDPDTaitwaterIsothermal::init_one (int i, int j) { + if (setflag[i][j] == 0) + error->all(FLERR,"Not all pair sph/taitwater/morris coeffs are not set"); + + cut[j][i] = cut[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- */ + +double PairSDPDTaitwaterIsothermal::single (int /*i*/, int /*j*/, int /*itype*/, + int /*jtype*/, double /*rsq*/, double /*factor_coul*/, + double /*factor_lj*/, double &fforce) { + fforce = 0.0; + + return 0.0; +} diff --git a/src/USER-SDPD/pair_sdpd_taitwater_isothermal.h b/src/USER-SDPD/pair_sdpd_taitwater_isothermal.h new file mode 100644 index 0000000000..afa21c4132 --- /dev/null +++ b/src/USER-SDPD/pair_sdpd_taitwater_isothermal.h @@ -0,0 +1,60 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(sdpd/taitwater/isothermal,PairSDPDTaitwaterIsothermal) + +#else + +#ifndef LMP_PAIR_SDPD_TAITWATER_MORRIS_ISOTHERMAL_H +#define LMP_PAIR_SDPD_TAITWATER_MORRIS_ISOTHERMAL_H + +#include "pair.h" +#ifdef USE_ZEST +#include +#include "zest.hpp" +#endif + +namespace LAMMPS_NS { + +class PairSDPDTaitwaterIsothermal : public Pair { + public: + PairSDPDTaitwaterIsothermal (class LAMMPS *); + virtual ~PairSDPDTaitwaterIsothermal (); + virtual void compute (int, int); + void settings (int, char **); + void coeff (int, char **); + virtual double init_one (int, int); + virtual double single (int, int, int, int, double, double, double, double &); + + protected: + double viscosity, temperature; + double *rho0, *soundspeed, *B; + double **cut; + + void allocate (); + + unsigned int seed; +#ifdef USE_ZEST + std::mt19937_64 generator; + Ziggurat gaussian; +#else + class RanMars *random; +#endif +}; + +} + +#endif +#endif