Merge branch 'master' into collected-small-fixes

find_path (NETCDF_INCLUDE_DIR netcdf.h

  HINTS "${NETCDF_DIR}/include")

mark_as_advanced (NETCDF_INCLUDE_DIR)

set (NETCDF_C_INCLUDE_DIRS ${NETCDF_INCLUDE_DIR})

string(REGEX REPLACE "/include/?$" ""

  NETCDF_LIB_HINT ${NETCDF_INCLUDE_DIR})

find_library (NETCDF_LIBRARY NAMES netcdf

  HINTS "${NETCDF_DIR}/lib")

  HINTS "${NETCDF_DIR}" "${NETCDF_LIB_HINT}" PATH_SUFFIXES lib lib64)

mark_as_advanced (NETCDF_LIBRARY)

set (NETCDF_C_LIBRARIES ${NETCDF_LIBRARY})

# source: https://ftp.space.dtu.dk/pub/Ioana/pism0.6.1-10/CMake/FindPNetCDF.cmake

# license: GPL v3 (https://ftp.space.dtu.dk/pub/Ioana/pism0.6.1-10/COPYING)

#

# - Find PNetCDF

# Find the native PNetCDF includes and library

#

#  PNETCDF_INCLUDES    - where to find netcdf.h, etc

#  PNETCDF_LIBRARIES   - Link these libraries when using NetCDF

#  PNETCDF_FOUND       - True if PNetCDF was found

#

# Normal usage would be:

#  find_package (PNetCDF REQUIRED)

#  target_link_libraries (uses_pnetcdf ${PNETCDF_LIBRARIES})

if (PNETCDF_INCLUDES AND PNETCDF_LIBRARIES)

  # Already in cache, be silent

  set (PNETCDF_FIND_QUIETLY TRUE)

endif (PNETCDF_INCLUDES AND PNETCDF_LIBRARIES)

find_path (PNETCDF_INCLUDES pnetcdf.h

  HINTS "${PNETCDF_ROOT}/include" "$ENV{PNETCDF_ROOT}/include")

string(REGEX REPLACE "/include/?$" ""

  PNETCDF_LIB_HINT ${PNETCDF_INCLUDES})

find_library (PNETCDF_LIBRARIES

  NAMES pnetcdf

  HINTS ${PNETCDF_LIB_HINT} PATH_SUFFIXES lib lib64)

if ((NOT PNETCDF_LIBRARIES) OR (NOT PNETCDF_INCLUDES))

  message(STATUS "Trying to find PNetCDF using LD_LIBRARY_PATH (we're desperate)...")

  file(TO_CMAKE_PATH "$ENV{LD_LIBRARY_PATH}" LD_LIBRARY_PATH)

  find_library(PNETCDF_LIBRARIES

    NAMES pnetcdf

    HINTS ${LD_LIBRARY_PATH})

  if (PNETCDF_LIBRARIES)

    get_filename_component(PNETCDF_LIB_DIR ${PNETCDF_LIBRARIES} PATH)

    string(REGEX REPLACE "/(lib|lib64)/?$" "/include"

      PNETCDF_H_HINT ${PNETCDF_LIB_DIR})

    find_path (PNETCDF_INCLUDES pnetcdf.h

      HINTS ${PNETCDF_H_HINT}

      DOC "Path to pnetcdf.h")

  endif()

endif()

# handle the QUIETLY and REQUIRED arguments and set PNETCDF_FOUND to TRUE if

# all listed variables are TRUE

include (FindPackageHandleStandardArgs)

find_package_handle_standard_args (PNetCDF DEFAULT_MSG PNETCDF_LIBRARIES PNETCDF_INCLUDES)

mark_as_advanced (PNETCDF_LIBRARIES PNETCDF_INCLUDES)

if(PKG_USER-NETCDF)

  find_package(NetCDF REQUIRED)

  # USER-NETCDF can use NetCDF, Parallel NetCDF (PNetCDF), or both. At least one necessary.

  # NetCDF library enables dump sytle "netcdf", while PNetCDF enables dump style "netcdf/mpiio"

  find_package(NetCDF)

  if(NETCDF_FOUND)

    find_package(PNetCDF)

  else(NETCDF_FOUND)

    find_package(PNetCDF REQUIRED)

  endif(NETCDF_FOUND)

  if(NETCDF_FOUND)

    include_directories(${NETCDF_INCLUDE_DIRS})

    list(APPEND LAMMPS_LINK_LIBS ${NETCDF_LIBRARIES})

  add_definitions(-DLMP_HAS_NETCDF -DNC_64BIT_DATA=0x0020)

    add_definitions(-DLMP_HAS_NETCDF)

  endif(NETCDF_FOUND)

  if(PNETCDF_FOUND)

    include_directories(${PNETCDF_INCLUDES})

    list(APPEND LAMMPS_LINK_LIBS ${PNETCDF_LIBRARIES})

    add_definitions(-DLMP_HAS_PNETCDF)

  endif(PNETCDF_FOUND)

  add_definitions(-DNC_64BIT_DATA=0x0020)

endif()

Here, *a* is the radius of the contact region, given by :math:`a =\sqrt{R\delta}`

for all normal contact models, except for *jkr*\ , where it is given

implicitly by :math:`\delta = a^2/R - 2\sqrt{\pi \gamma a/E}`, see

discussion above. To match the Mindlin solution, one should set :math:`k_t = 8G`, where :math:`G` is the shear modulus, related to Young's modulus

discussion above. To match the Mindlin solution, one should set :math:`k_t = 4G/(2-\nu)`, where :math:`G` is the shear modulus, related to Young's modulus

:math:`E` by :math:`G = E/(2(1+\nu))`, where :math:`\nu` is Poisson's ratio. This

can also be achieved by specifying *NULL* for :math:`k_t`, in which case a

normal contact model that specifies material parameters :math:`E` and

#Pour composite granular particles on flat wall

newton		on

atom_style	sphere

atom_modify	map array sort 0 0

thermo_modify	flush yes

units 		si  

variable	minrad equal 0.5

variable	maxrad equal 1.4

variable 	skin equal 0.3*${maxrad}

boundary	p p f

region		reg block 0 20 0 20 0 200 units box

create_box	1 reg

fix             prop all property/atom mol ghost yes

variable	dumpfreq equal 1000

variable	logfreq equal 1000

pair_style	gran/hooke/history 4e5 NULL 1e2 NULL 0.5 0

pair_coeff	* *

timestep 	0.0001

group		particles type 1

atom_modify     first particles

neighbor        ${skin} bin

group		rigid type 1

neigh_modify	every 1 delay 0 check yes exclude molecule/intra all 

thermo 		${logfreq}

thermo_style	custom step cpu atoms ke

thermo_modify	flush yes lost warn

comm_modify	vel yes cutoff 3

molecule	mymol molecule.data

region		pourreg block 5 15 5 15 80 100 side in units box

#Note: in versions prior to 1/2020, the 'disable' keyword to fix/gravity

# and the 'gravity' keyword to fix rigid/small were not available.

# These settings produce undesirable behavior, where gravity can induce

# torque on rigid bodies.

#fix		gravfix all gravity 9.8 vector 0 0 -1 #disable

#fix		rigidfix all rigid/small molecule mol mymol #gravity gravfix

#The correct behavior is recovered with the following settings:

fix		gravfix all gravity 9.8 vector 0 0 -1 disable

fix		rigidfix all rigid/small molecule mol mymol gravity gravfix

fix		pourfix all pour 5 0 1234 region pourreg mol mymol rigid rigidfix

fix		zwall all wall/gran hooke/history 4000.0 NULL 100.0 NULL 0.5 0 zplane 0.1 NULL

#dump		1 all custom 1000 molecule_pour.dump id type mass radius x y z fx fy fz

run		100000