diff --git a/examples/KAPPA/README b/examples/KAPPA/README index ee70766d9..d025797bf 100644 --- a/examples/KAPPA/README +++ b/examples/KAPPA/README @@ -1,104 +1,104 @@ This directory has 5 scripts that compute the thermal conductivity (kappa) of a Lennard-Jones fluid using 5 different methods. See the discussion in Section 6.20 of the manual for an overview of the methods and pointers to doc pages for the commands which implement them. Citations for the various methods can also be found in the manaul. These scripts are provided for illustration purposes. No guarantee is made that the systems are fully equilibrated or that the runs are long enough to generate good statistics and highly accurate results. These scripts could easily be adapted to work with solids as well. ------------- These are the 5 methods for computing thermal conductivity. The first 4 are non-equilibrium methods; the last is an equilibrium method. in.langevin = thermostat 2 regions at different temperatures via fix langevin in.heat = add/subtract energy to 2 regions via fix heat in.ehex = add/subtract energy to 2 regions via fix ehex in.mp = use fix thermal/conductivity and the Muller-Plathe method in.heatflux = use compute heat/flux and the Green-Kubo method The NEMD systems have 8000 atoms with a box length 2x larger in z, the non-equilibrium direction. The G-K system has 4000 atoms and a cubic box; it also needs to be run longer to generate good statistics. The scripts were all run on 8 processors. They all run in a minute or so and produce the accompanying log files and profile files (for temperature or heat-flux). The state point of the LJ fluid is rho* = 0.6, T* = 1.35, and Rcut = 2.5 sigma. This was chosen to agree with a 1986 paper by D Evans in Phys Rev A, 34, p 1449, where he computed the thermal conductivity of a small 108-atom system using a thermostatting method. Fig 1 in the paper shows his simulations produced a kappa of around 3.4 for this system, in agreement with an experimental data point as well. ------------- Here is how to extract Kappa from the log file output for each method. The NEMD methods use the formula kappa = dQ * dZ/dTemp where dQ = energy flux, and dTemp/dZ = temperature gradient. (1) in.langevin dQ = 8000 * 0.5*(0.905+0.947) / 100 / 18.82^2 / 2 8000 atoms 0.5*(0.905+0.947) = from log file = ave of total in/out energy for 2 regions normalized by # of atoms 100 = 20,000 steps at 0.005 tau timestep = run time in tau xy box area = 18.82^2 divide by 2 since energy flux goes in 2 directions due to periodic z dTemp = 0.578 from log file for average Temp difference between 2 regions dZ = 18.82 Kappa = 3.41 (2) in.heat dQ = (100*100) / 100 / 18.82^2 / 2 100*100 = 100 (time in tau) * 100 (energy delta specified in fix heat) 100 = 20,000 steps at 0.005 tau timestep = run time in tau xy box area = 18.82^2 divide by 2 since energy flux goes in 2 directions due to periodic z dTemp = 0.748 from log file for average Temp difference between 2 regions dZ = 18.82 Kappa = 3.55 (3) in.ehex dQ = (100*100) / 100 / 18.82^2 / 2 100*100 = 100 (time in tau) * 100 (energy delta specified in fix heat) 100 = 20,000 steps at 0.005 tau timestep = run time in tau xy box area = 18.82^2 divide by 2 since energy flux goes in 2 directions due to periodic z dTemp = 0.770 from log file for average Temp difference between 2 regions dZ = 18.82 Kappa = 3.45 (4) in.mp dQ = 15087 / 100 / 18.82^2 / 2 - 15087 = cummulative delta energy, tallied by fix thermal/conductivity + 15087 = cumulative delta energy, tallied by fix thermal/conductivity 100 = 20,000 steps at 0.005 tau timestep = run time in tau xy box area = 18.82^2 divide by 2 since energy flux goes in 2 directions due to periodic z dTemp = 1.16 from log file for average Temp difference between 2 regions dZ = 18.82 Kappa = 3.45 (5) in.heatflux kappa is computed directly within the script, by performing a time integration of the formulas discussed on the compute heat/flux doc page - the resulting value prints at the end of the run and is in the log file Kappa = 3.78 diff --git a/examples/USER/lb/confined_colloid/in.confined_colloids b/examples/USER/lb/confined_colloid/in.confined_colloids index 37cd52e62..b7b8958e9 100755 --- a/examples/USER/lb/confined_colloid/in.confined_colloids +++ b/examples/USER/lb/confined_colloid/in.confined_colloids @@ -1,108 +1,108 @@ #===========================================================================# # System of colloidal particles under shear. # # # # Run consists of 10x12x4 particles, each composed of 3613 particle nodes # # (3612 particles forming a spherical shell, and 1 central particle). # # 280 x 280 x 101 lattice-Boltzmann grid sites. # # # # This simulation is used to illustrate the simulation time of a realistic # # implementation of the lb_fluid fix. # # The data file "confinedcolloids.dat" is quite large and so is not # # included here. It can be obtained from: # # http://www.apmaths.uwo.ca/~cdennist/confinedcolloids.dat.gz # # # # Sample output from this run can be found in the file: # # results64.out # #===========================================================================# units micro dimension 3 boundary p p f atom_style molecular read_data confinedcolloids.dat mass * 0.00010287 #---------------------------------------------------------------------------- # Need a neighbor bin size smaller than the lattice-Boltzmann grid spacing # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. # The arguments for neigh_modify have been set to "delay 0 every 1", again # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. However, these values can likely # be somewhat increased without issue. If a problem does arise (a particle # is outside of its processors LB grid) an error message is printed and # the simulation is terminated. #---------------------------------------------------------------------------- neighbor 0.03 bin neigh_modify delay 0 every 1 neigh_modify exclude type 2 2 neigh_modify exclude type 2 1 #---------------------------------------------------------------------------- # ForceAtoms are the particles at the center of each colloidal object which # do not interact with the fluid, but are used to implement the hard-sphere # interactions. #---------------------------------------------------------------------------- group ForceAtoms type 1 #---------------------------------------------------------------------------- # FluidAtoms are the particles representing the surface of the colloidal # object which do interact with the fluid. #---------------------------------------------------------------------------- group FluidAtoms type 2 #---------------------------------------------------------------------------- # Implement a hard-sphere interaction between the particles at the center of # each colloidal object (use a truncated and shifted Lennard-Jones # potential). #---------------------------------------------------------------------------- pair_style lj/cut 1.572 pair_coeff * * 0.0 0.0 1.572 pair_coeff 1 1 10.0 1.400492785 1.572 pair_modify shift yes timestep 0.0006 #--------------------------------------------------------------------------- # Create a lattice-Boltzmann fluid covering the simulation domain. # This fluid feels a force due to the particles specified through FluidAtoms -# (however, this fix does not explicity apply a force back on to these +# (however, this fix does not explicitly apply a force back on to these # particles...this is accomplished through the use of the viscous_lb fix). # Use the standard LB integration scheme, a fluid density = 1.0, # fluid viscosity = 1.0, lattice spacing dx=0.06, and mass unit, dm=0.00003. # Use the default method to calculate the interaction force between the # particles and the fluid. This calculation makes use of the surface area # of the composite object represented by each particle node. Since this is # not equal to dx*dx, it is input through the setArea keyword (i.e. # particles of type 2 correspond to a surface area of 0.0018337299). # Implement walls moving at speeds of 20.0 in opposite directions. #---------------------------------------------------------------------------- fix 1 FluidAtoms lb/fluid 1 1 1.0 1.0 dx 0.06 dm 0.00003 setArea 2 0.0018337299 zwall_velocity -20.0 20.0 #---------------------------------------------------------------------------- # Apply the force due to the fluid onto the FluidAtoms particles (again, # these atoms represent the surface of the colloidal object, which should # interact with the fluid). #---------------------------------------------------------------------------- fix 2 FluidAtoms lb/viscous #---------------------------------------------------------------------------- # Each colloidal object (spherical shell of particles and central particle) # is specified as a separate molecule in the confinedcolloids.dat data # file. Integrate the motion of each of these sets of particles as rigid # objects which move and rotate together. #---------------------------------------------------------------------------- fix 3 all rigid molecule #---------------------------------------------------------------------------- # Implement a repulsive interaction between the ForceAtoms particles, and the # upper and lower z-walls. (A truncated and shifted Lennard-Jones potential # is used). #---------------------------------------------------------------------------- fix wallhi ForceAtoms wall/lj126 zhi 5.88 20.0 0.8071542386 0.906 units box fix walllo ForceAtoms wall/lj126 zlo 0.0 20.0 0.8071542386 0.906 units box #dump ParticleTracking ForceAtoms custom 50 test.track id mol x y z vx vy vz run 400 diff --git a/examples/USER/lb/microrheology/in.microrheology_default_gamma b/examples/USER/lb/microrheology/in.microrheology_default_gamma index 4937a2cb2..8c3b684ff 100755 --- a/examples/USER/lb/microrheology/in.microrheology_default_gamma +++ b/examples/USER/lb/microrheology/in.microrheology_default_gamma @@ -1,119 +1,119 @@ #===========================================================================# # 2 particle microrheology test # # # # Run consists of 2 colloidal particles undergoing Brownian motion in a # # thermal lattice-Boltzmann fluid. # # # # Here, gamma (used in the calculation of the particle-fluid interaction # # force) is calculated by default. Thus, the colloidal objects will have # # a slightly larger "hydrodynamic" radii than given by the placement of # # the particle nodes. # # # # Sample output from this run can be found in the file: # # 'microrheology_setgamma.out' # #===========================================================================# units nano dimension 3 boundary p p p atom_style molecular read_data data.two #---------------------------------------------------------------------------- # Need a neighbor bin size smaller than the lattice-Boltzmann grid spacing # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. # The arguments for neigh_modify have been set to "delay 0 every 1", again # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. However, these values can likely # be somewhat increased without issue. If a problem does arise (a particle # is outside of its processors LB grid) an error message is printed and # the simulation is terminated. #---------------------------------------------------------------------------- neighbor 0.3 bin neigh_modify delay 0 every 1 neigh_modify exclude type 2 2 neigh_modify exclude type 2 1 #---------------------------------------------------------------------------- # Implement a hard-sphere interaction between the particles at the center of # each colloidal object (use a truncated and shifted Lennard-Jones # potential). #---------------------------------------------------------------------------- pair_style lj/cut 5.88 pair_coeff * * 0.0 0.0 5.88 pair_coeff 1 1 100.0 5.238484463 5.88 pair_modify shift yes mass * 0.0002398 timestep 0.00025 #---------------------------------------------------------------------------- # ForceAtoms are the particles at the center of each colloidal object which # do not interact with the fluid, but are used to implement the hard-sphere # interactions. # FluidAtoms are the particles representing the surface of the colloidal # object which do interact with the fluid. #---------------------------------------------------------------------------- group ForceAtoms type 1 group FluidAtoms type 2 #--------------------------------------------------------------------------- # Create a lattice-Boltzmann fluid covering the simulation domain. # This fluid feels a force due to the particles specified through FluidAtoms -# (however, this fix does not explicity apply a force back on to these +# (however, this fix does not explicitly apply a force back on to these # particles...this is accomplished through the use of the viscous_lb fix). # Use the standard LB integration scheme, a fluid viscosity = 1.0, fluid # density= 0.0009982071, lattice spacing dx=1.2, and mass unit, dm=0.003. # Use the default method to calculate the interaction force between the # particles and the fluid. This calculation requires the surface area # of the composite object represented by each particle node. By default # this area is assumed equal to dx*dx; however, since this is not the case # here, it is input through the setArea keyword (i.e. particles of type 2 # correspond to a surface area of 0.3015928947). # Use the trilinear interpolation stencil to distribute the force from # a given particle onto the fluid mesh. # Use a thermal lattice-Boltzmann fluid (temperature 300K, random number # seed=2762). This enables the particles to undergo Brownian motion in # the fluid. #---------------------------------------------------------------------------- fix 1 FluidAtoms lb/fluid 1 1 1.0 0.0009982071 setArea 2 0.3015928947 dx 1.2 dm 0.003 trilinear noise 300.0 2762 #---------------------------------------------------------------------------- # Apply the force due to the fluid onto the FluidAtoms particles (again, # these atoms represent the surface of the colloidal object, which should # interact with the fluid). #---------------------------------------------------------------------------- fix 2 FluidAtoms lb/viscous #---------------------------------------------------------------------------- # Each colloidal object (spherical shell of particles and central particle) # is specified as a separate molecule in the confinedcolloids.dat data # file. Integrate the motion of these sets of particles as rigid objects # which each move and rotate together. #---------------------------------------------------------------------------- fix 3 all rigid molecule #---------------------------------------------------------------------------- # To ensure that numerical errors do not lead to a buildup of momentum in the # system, the momentum_lb fix is used every 10000 timesteps to zero out the # total (particle plus fluid) momentum in the system. #---------------------------------------------------------------------------- fix 4 all lb/momentum 10000 linear 1 1 1 #---------------------------------------------------------------------------- # Create variables containing the positions of the central atoms (these # values should correspond to the center of mass of each composite # colloidal particle), and output these quantities to the screen. #---------------------------------------------------------------------------- variable x1 equal x[1] variable y1 equal y[1] variable z1 equal z[1] variable x2 equal x[242] variable y2 equal y[242] variable z2 equal z[242] thermo_style custom v_x1 v_y1 v_z1 v_x2 v_y2 v_z2 thermo 1 run 2000000000 diff --git a/examples/USER/lb/microrheology/in.microrheology_set_gamma b/examples/USER/lb/microrheology/in.microrheology_set_gamma index 7b84bdaf6..1f744220f 100755 --- a/examples/USER/lb/microrheology/in.microrheology_set_gamma +++ b/examples/USER/lb/microrheology/in.microrheology_set_gamma @@ -1,113 +1,113 @@ #===========================================================================# # 2 particle microrheology test # # # # Run consists of 2 colloidal particles undergoing Brownian motion in a # # thermal lattice-Boltzmann fluid. # # # # Here, gamma (used in the calculation of the particle-fluid interaction # # force) is set by the user (gamma = 1.4692 for this simulation...this # # value has been calibrated a priori through simulations of the drag # # force acting on a single particle of the same radius). # # # # Sample output from this run can be found in the file: # # 'microrheology_setgamma.out' # #===========================================================================# units nano dimension 3 boundary p p p atom_style molecular read_data data.two #---------------------------------------------------------------------------- # Need a neighbor bin size smaller than the lattice-Boltzmann grid spacing # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. # The arguments for neigh_modify have been set to "delay 0 every 1", again # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. However, these values can likely # be somewhat increased without issue. If a problem does arise (a particle # is outside of its processors LB grid) an error message is printed and # the simulation is terminated. #---------------------------------------------------------------------------- neighbor 0.3 bin neigh_modify delay 0 every 1 neigh_modify exclude type 2 2 neigh_modify exclude type 2 1 #---------------------------------------------------------------------------- # Implement a hard-sphere interaction between the particles at the center of # each colloidal object (use a truncated and shifted Lennard-Jones # potential). #---------------------------------------------------------------------------- pair_style lj/cut 5.88 pair_coeff * * 0.0 0.0 5.88 pair_coeff 1 1 100.0 5.238484463 5.88 pair_modify shift yes mass * 0.0002398 timestep 0.00045 #---------------------------------------------------------------------------- # ForceAtoms are the particles at the center of each colloidal object which # do not interact with the fluid, but are used to implement the hard-sphere # interactions. # FluidAtoms are the particles representing the surface of the colloidal # object which do interact with the fluid. #---------------------------------------------------------------------------- group ForceAtoms type 1 group FluidAtoms type 2 #--------------------------------------------------------------------------- # Create a lattice-Boltzmann fluid covering the simulation domain. # This fluid feels a force due to the particles specified through FluidAtoms -# (however, this fix does not explicity apply a force back on to these +# (however, this fix does not explicitly apply a force back on to these # particles...this is accomplished through the use of the rigid_pc_sphere # fix). # Use the LB integration scheme of Ollila et. al. (for stability reasons, # this integration scheme should be used when a large user set value for # gamma is specified), a fluid viscosity = 1.0, fluid density= 0.0009982071, # value for gamma=1.4692, lattice spacing dx=1.2, and mass unit, dm=0.003. # Use a thermal lattice-Boltzmann fluid (temperature 300K, random number # seed=2762). This enables the particles to undergo Brownian motion in # the fluid. #---------------------------------------------------------------------------- fix 1 FluidAtoms lb/fluid 1 2 1.0 0.0009982071 setGamma 1.4692 dx 1.2 dm 0.003 noise 300.0 2762 #---------------------------------------------------------------------------- # Apply the force from the fluid to the particles, and integrate their # motion, constraining them to move and rotate together as a single rigid # spherical object. # Since both the ForceAtoms (central atoms), and the FluidAtoms (spherical # shell) should move and rotate together, this fix is applied to all of # the atoms in the system. However, since the central atoms should not # feel a force due to the fluid, they are excluded from the fluid force # calculation through the use of the 'innerNodes' keyword. # NOTE: This fix should only be used when the user specifies a value for # gamma (through the setGamma keyword) in the lb_fluid fix. #---------------------------------------------------------------------------- fix 2 all lb/rigid/pc/sphere molecule innerNodes ForceAtoms #---------------------------------------------------------------------------- # To ensure that numerical errors do not lead to a buildup of momentum in the # system, the momentum_lb fix is used every 10000 timesteps to zero out the # total (particle plus fluid) momentum in the system. #---------------------------------------------------------------------------- fix 3 all lb/momentum 10000 linear 1 1 1 #---------------------------------------------------------------------------- # Create variables containing the positions of the central atoms (these # values should correspond to the center of mass of each composite # colloidal particle), and output these quantities to the screen. #---------------------------------------------------------------------------- variable x1 equal x[1] variable y1 equal y[1] variable z1 equal z[1] variable x2 equal x[242] variable y2 equal y[242] variable z2 equal z[242] thermo_style custom v_x1 v_y1 v_z1 v_x2 v_y2 v_z2 thermo 1 run 2000000000 diff --git a/examples/USER/lb/planewall/in.planewall_default_gamma b/examples/USER/lb/planewall/in.planewall_default_gamma index bf0cfcb3c..83f42ebdd 100755 --- a/examples/USER/lb/planewall/in.planewall_default_gamma +++ b/examples/USER/lb/planewall/in.planewall_default_gamma @@ -1,99 +1,99 @@ #===========================================================================# # Rigid sphere freely moving near a stationary plane wall in a system # # undergoing shear flow. # # Every 10 time steps the center of mass velocity and angular velocity of # # the sphere are printed to the screen. # # # # Here, gamma (used in the calculation of the particle-fluid interaction # # force) is calculated by default. Thus, the colloidal objects will have # # a slightly larger "hydrodynamic" radii than given by the placement of # # the particle nodes. # # # # Sample output from this run can be found in the file: # # 'wall_defaultgamma.out' # #===========================================================================# units micro dimension 3 boundary p p f atom_style atomic #---------------------------------------------------------------------------- # Need a neighbor bin size smaller than the lattice-Boltzmann grid spacing # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. # The arguments for neigh_modify have been set to "delay 0 every 1", again # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. However, these values can likely # be somewhat increased without issue. If a problem does arise (a particle # is outside of its processors LB grid) an error message is printed and # the simulation is terminated. #---------------------------------------------------------------------------- neighbor 1.0 bin neigh_modify delay 0 every 1 read_data data.one_radius16d2 #---------------------------------------------------------------------------- # None of the particles interact with one another. #---------------------------------------------------------------------------- pair_style lj/cut 2.45 pair_coeff * * 0.0 0.0 2.45 neigh_modify exclude type 1 1 mass * 100.0 timestep 3.0 group sphere1 id <> 1 320 #---------------------------------------------------------------------------- # Colloidal particle is initially stationary. #---------------------------------------------------------------------------- velocity all set 0.0 0.0 0.0 units box #---------------------------------------------------------------------------- # Create a lattice-Boltzmann fluid covering the simulation domain. # All of the particles in the simulation apply a force to the fluid. -# (however, this fix does not explicity apply a force back on to these +# (however, this fix does not explicitly apply a force back on to these # particles...this is accomplished through the use of the viscous_lb fix. # Use the standard LB integration scheme, a fluid density = 1.0, # fluid viscosity = 1.0, lattice spacing dx=4.0, and mass unit, dm=10.0. # Use the default method to calculate the interaction force between the # particles and the fluid. This calculation requires the surface area # of the composite object represented by each particle node. By default # this area is assumed equal to dx*dx; however, since this is not the case # here, it is input through the setArea keyword (i.e. particles of type 1 # correspond to a surface area of 10.3059947). # Use the trilinear interpolation stencil to distribute the force from # a given particle onto the fluid mesh. # Create shear in the system, by giving the upper z-wall a velocity of 0.0001 # along the y-direction, while keeping the lower z-wall stationary. #----------------------------------------------------------------------------- fix 1 all lb/fluid 1 1 1.0 1.0 setArea 1 10.3059947 dx 4.0 dm 10.0 trilinear zwall_velocity 0.0 0.0001 #---------------------------------------------------------------------------- # Apply the force due to the fluid onto the particles. #---------------------------------------------------------------------------- fix 2 all lb/viscous #---------------------------------------------------------------------------- # Integrate the motion of the particles, constraining them to move and # rotate together as a single rigid spherical object. #---------------------------------------------------------------------------- fix 3 all rigid group 1 sphere1 #---------------------------------------------------------------------------- # Create variables for the center-of-mass and angular velocities, and output # these quantities to the screen. #---------------------------------------------------------------------------- variable vx equal vcm(all,x) variable vy equal vcm(all,y) variable vz equal vcm(all,z) variable omegax equal omega(all,x) variable omegay equal omega(all,y) variable omegaz equal omega(all,z) thermo_style custom v_vx v_vy v_vz v_omegax v_omegay v_omegaz thermo 10 run 200000 diff --git a/examples/USER/lb/planewall/in.planewall_set_gamma b/examples/USER/lb/planewall/in.planewall_set_gamma index a048411b0..47d8266a1 100755 --- a/examples/USER/lb/planewall/in.planewall_set_gamma +++ b/examples/USER/lb/planewall/in.planewall_set_gamma @@ -1,92 +1,92 @@ #===========================================================================# # Rigid sphere freely moving near a stationary plane wall in a system # # undergoing shear flow. # # Every 10 time steps the center of mass velocity and angular velocity of # # the sphere are printed to the screen. # # # # Here, gamma (used in the calculation of the particle-fluid interaction # # force) is set by the user (gamma = 13.655 for this simulation...this # # value has been calibrated a priori through simulations of the drag # # force acting on a single particle of the same radius). # # # # Sample output from this run can be found in the file: # # 'wall_setgamma.out' # #===========================================================================# units micro dimension 3 boundary p p f atom_style atomic #---------------------------------------------------------------------------- # Need a neighbor bin size smaller than the lattice-Boltzmann grid spacing # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. # The arguments for neigh_modify have been set to "delay 0 every 1", again # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. However, these values can likely # be somewhat increased without issue. If a problem does arise (a particle # is outside of its processors LB grid) an error message is printed and # the simulation is terminated. #---------------------------------------------------------------------------- neighbor 1.0 bin neigh_modify delay 0 every 1 read_data data.one_radius16d2 #---------------------------------------------------------------------------- # None of the particles interact with one another. #---------------------------------------------------------------------------- pair_style lj/cut 2.45 pair_coeff * * 0.0 0.0 2.45 neigh_modify exclude type 1 1 mass * 100.0 timestep 4.0 group sphere1 id <> 1 320 #---------------------------------------------------------------------------- # Colloidal particle is initially stationary. #---------------------------------------------------------------------------- velocity all set 0.0 0.0 0.0 units box #---------------------------------------------------------------------------- # Create a lattice-Boltzmann fluid covering the simulation domain. # All of the particles in the simulation apply a force to the fluid. -# (however, this fix does not explicity apply a force back on to these +# (however, this fix does not explicitly apply a force back on to these # particles...this is accomplished through the use of the rigid_pc_sphere # fix). # Use the LB integration scheme of Ollila et. al. (for stability reasons, # this integration scheme should be used when a large user set value for # gamma is specified), a fluid density = 1.0, fluid viscosity = 1.0, value # for gamma=13.655, lattice spacing dx=4.0, and mass unit, dm=10.0. # Create shear in the system, by giving the upper z-wall a velocity of 0.0001 # along the y-direction, while keeping the lower z-wall stationary. #----------------------------------------------------------------------------- fix 1 all lb/fluid 1 2 1.0 1.0 setGamma 13.655 dx 4.0 dm 10.0 zwall_velocity 0.0 0.0001 #---------------------------------------------------------------------------- # Apply the force from the fluid to the particles, and integrate their # motion, constraining them to move and rotate together as a single rigid # spherical object. # NOTE: This fix should only be used when the user specifies a value for # gamma (through the setGamma keyword) in the lb_fluid fix. #---------------------------------------------------------------------------- fix 2 all lb/rigid/pc/sphere group 1 sphere1 #---------------------------------------------------------------------------- # Create variables for the center-of-mass and angular velocities, and output # these quantities to the screen. #---------------------------------------------------------------------------- variable vx equal vcm(all,x) variable vy equal vcm(all,y) variable vz equal vcm(all,z) variable omegax equal omega(all,x) variable omegay equal omega(all,y) variable omegaz equal omega(all,z) thermo_style custom v_vx v_vy v_vz v_omegax v_omegay v_omegaz thermo 10 run 200000 diff --git a/examples/USER/lb/polymer/in.polymer_default_gamma b/examples/USER/lb/polymer/in.polymer_default_gamma index 4344dc0ef..cd48a8ccf 100755 --- a/examples/USER/lb/polymer/in.polymer_default_gamma +++ b/examples/USER/lb/polymer/in.polymer_default_gamma @@ -1,107 +1,107 @@ #===========================================================================# # polymer test # # # # Run consists of a lone 32-bead coarse-grained polymer # # undergoing Brownian motion in thermal lattice-Boltzmann fluid. # # # # Here, gamma (used in the calculation of the monomer-fluid interaction # # force) is set by the user (gamma = 0.03 for this simulation...this # # value has been calibrated a priori through simulations of the drag # # force acting on a single particle of the same radius). # # Sample output from this run can be found in the file: # # 'dump.polymer.lammpstrj' # # and viewed using, e.g., the VMD software. # # # #===========================================================================# units nano dimension 3 boundary p p p atom_style hybrid molecular special_bonds fene read_data data.polymer #---------------------------------------------------------------------------- # Need a neighbor bin size smaller than the lattice-Boltzmann grid spacing # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. #---------------------------------------------------------------------------- neighbor 0.5 bin neigh_modify delay 0 every 1 check yes neigh_modify exclude type 2 2 neigh_modify exclude type 2 1 #---------------------------------------------------------------------------- # Implement a hard-sphere interaction between the particles at the center of # each monomer (use a truncated and shifted Lennard-Jones potential). #---------------------------------------------------------------------------- bond_style fene bond_coeff 1 60.0 2.25 4.14195 1.5 pair_style lj/cut 1.68369 pair_coeff 1 1 4.14195 1.5 1.68369 pair_coeff 1 2 4.14195 1.5 1.68369 pair_coeff 2 2 0 1.0 mass * 0.000000771064 timestep 0.00003 #---------------------------------------------------------------------------- # ForceAtoms are the particles at the center of each monomer which # do not interact with the fluid, but are used to implement the hard-sphere # interactions. # FluidAtoms are the particles representing the surface of the monomer # which do interact with the fluid. Monomer surface is shell of radius 0.7 #---------------------------------------------------------------------------- group ForceAtoms type 1 group FluidAtoms type 2 #--------------------------------------------------------------------------- # Create a lattice-Boltzmann fluid covering the simulation domain. # This fluid feels a force due to the particles specified through FluidAtoms -# (however, this fix does not explicity apply a force back on to these +# (however, this fix does not explicitly apply a force back on to these # particles. This is accomplished through the use of the lb/viscous # fix). # Uses the standard LB integration scheme, fluid viscosity = 0.023333333, # fluid density= 0.0000166368, lattice spacing dx=1.0, and mass unit, # dm=0.0000166368. # Use the default method to calculate the interaction force between the # particles and the fluid. This calculation requires the surface area # of the composite object represented by each particle node. By default # this area is assumed equal to dx*dx; however, since this is not the case # here, it is input through the setArea keyword (i.e. particles of type 2 # correspond to a surface area of 0.2025=4 Pi R^2/N ). # Use the trilinear interpolation stencil to distribute the force from # a given particle onto the fluid mesh (results in a smaller hydrodynamic # radius than if the Peskin stencil is used). # Use a thermal lattice-Boltzmann fluid (temperature 300K, random number # seed=15003). This enables the particles to undergo Brownian motion in # the fluid. #---------------------------------------------------------------------------- fix 1 FluidAtoms lb/fluid 3 1 0.023333333 0.0000166368 setArea 2 0.20525 dx 1.0 dm 0.0000166368 noise 300.0 15003 #---------------------------------------------------------------------------- # Apply the force from the fluid to the particles, and integrate their # motion, constraining them to move and rotate together as a single rigid # spherical object. # Since both the ForceAtoms (central atoms), and the FluidAtoms (spherical # shell) should move and rotate together, this fix is applied to all of # the atoms in the system. However, since the central atoms should not # feel a force due to the fluid, they are excluded from the fluid force # calculation. #---------------------------------------------------------------------------- fix 2 FluidAtoms lb/viscous fix 3 all rigid molecule #---------------------------------------------------------------------------- # To ensure that numerical errors do not lead to a buildup of momentum in the # system, the momentum_lb fix is used every 10000 timesteps to zero out the # total (particle plus fluid) momentum in the system. #---------------------------------------------------------------------------- fix 4 all lb/momentum 10000 linear 1 1 1 #---------------------------------------------------------------------------- # Write position and velocity coordinates into a file every 2000 time steps. #---------------------------------------------------------------------------- dump 1 ForceAtoms custom 2000 dump.polymer_default_gamma.lammpstrj id x y z vx vy vz run 2000001 diff --git a/examples/USER/lb/polymer/in.polymer_setgamma b/examples/USER/lb/polymer/in.polymer_setgamma index b5108a60f..fd5aeaa83 100755 --- a/examples/USER/lb/polymer/in.polymer_setgamma +++ b/examples/USER/lb/polymer/in.polymer_setgamma @@ -1,105 +1,105 @@ #===========================================================================# # polymer test # # # # Run consists of a lone 32-bead coarse-grained polymer # # undergoing Brownian motion in thermal lattice-Boltzmann fluid. # # # # Here, gamma (used in the calculation of the monomer-fluid interaction # # force) is set by the user (gamma = 0.03 for this simulation...this # # value has been calibrated a priori through simulations of the drag # # force acting on a single particle of the same radius). # # Sample output from this run can be found in the file: # # 'dump.polymer.lammpstrj' # # and viewed using, e.g., the VMD software. # # # # Santtu Ollila # # santtu.ollila@aalto.fi # # Aalto University # # August 14, 2013 # #===========================================================================# units nano dimension 3 boundary p p p atom_style hybrid molecular special_bonds fene read_data data.polymer #---------------------------------------------------------------------------- # Need a neighbor bin size smaller than the lattice-Boltzmann grid spacing # to ensure that the particles belonging to a given processor remain inside # that processors lattice-Boltzmann grid. #---------------------------------------------------------------------------- neighbor 0.5 bin neigh_modify delay 0 every 1 check yes neigh_modify exclude type 2 2 neigh_modify exclude type 2 1 #---------------------------------------------------------------------------- # Implement a hard-sphere interaction between the particles at the center of # each monomer (use a truncated and shifted Lennard-Jones potential). #---------------------------------------------------------------------------- bond_style fene bond_coeff 1 60.0 2.25 4.14195 1.5 pair_style lj/cut 1.68369 pair_coeff 1 1 4.14195 1.5 1.68369 pair_coeff 1 2 4.14195 1.5 1.68369 pair_coeff 2 2 0 1.0 mass * 0.000000771064 timestep 0.00003 #---------------------------------------------------------------------------- # ForceAtoms are the particles at the center of each monomer which # do not interact with the fluid, but are used to implement the hard-sphere # interactions. # FluidAtoms are the particles representing the surface of the monomer # which do interact with the fluid. #---------------------------------------------------------------------------- group ForceAtoms type 1 group FluidAtoms type 2 #--------------------------------------------------------------------------- # Create a lattice-Boltzmann fluid covering the simulation domain. # This fluid feels a force due to the particles specified through FluidAtoms -# (however, this fix does not explicity apply a force back on to these +# (however, this fix does not explicitly apply a force back on to these # particles. This is accomplished through the use of the rigid_pc_sphere # fix). # Use the LB integration scheme of Ollila et. al. (for stability reasons, # this integration scheme should be used when a large user set value for # gamma is specified), a fluid viscosity = 0.023333333, # fluid density= 0.0000166368, # value for gamma=0.03, lattice spacing dx=1.0, and mass unit, dm=0.0000166368. # Use a thermal lattice-Boltzmann fluid (temperature 300K, random number # seed=15003). This enables the particles to undergo Brownian motion in # the fluid. #---------------------------------------------------------------------------- fix 1 FluidAtoms lb/fluid 5 1 0.023333333 0.0000166368 setGamma 0.03 dx 1.0 dm 0.0000166368 noise 300.0 15003 #---------------------------------------------------------------------------- # Apply the force from the fluid to the particles, and integrate their # motion, constraining them to move and rotate together as a single rigid # spherical object. # Since both the ForceAtoms (central atoms), and the FluidAtoms (spherical # shell) should move and rotate together, this fix is applied to all of # the atoms in the system. However, since the central atoms should not # feel a force due to the fluid, they are excluded from the force # calculation through the use of the 'innerNodes' keyword. # NOTE: This fix should only be used when the user specifies a value for # gamma (through the setGamma keyword) in the lb_fluid fix. #---------------------------------------------------------------------------- fix 2 all lb/rigid/pc/sphere molecule innerNodes ForceAtoms #---------------------------------------------------------------------------- # To ensure that numerical errors do not lead to a buildup of momentum in the # system, the momentum_lb fix is used every 10000 timesteps to zero out the # total (particle plus fluid) momentum in the system. #---------------------------------------------------------------------------- fix 3 all lb/momentum 10000 linear 1 1 1 #---------------------------------------------------------------------------- # Write position and velocity coordinates into a file every 2000 time steps. #---------------------------------------------------------------------------- dump 1 ForceAtoms custom 2000 dump.polymer_setgamma.lammpstrj id x y z vx vy vz run 2000001 diff --git a/examples/VISCOSITY/README b/examples/VISCOSITY/README index 88b874cd6..98ca9fbe3 100644 --- a/examples/VISCOSITY/README +++ b/examples/VISCOSITY/README @@ -1,91 +1,91 @@ This directory has 5 scripts that compute the viscosity (eta) of a Lennard-Jones fluid using 5 different methods. See the discussion in Section 6.21 of the manual for an overview of the methods and pointers to doc pages for the commands which implement them. Citations for the various methods can also be found in the manual. These scripts are provided for illustration purposes. No guarantee is made that the systems are fully equilibrated or that the runs are long enough to generate good statistics and highly accurate results. ------------- These are the 5 methods for computing viscosity. The first 3 are non-equilibrium methods; the last 2 are equilibrium methods. in.wall = move a wall to shear the fluid between two walls in.nemd = use fix deform and fix nvt/sllod to perform a NEMD shear simulation in.mp = use fix viscosity and the Muller-Plathe method in.gk = use the Green-Kubo method in.einstein = use the Einstein version of Green-Kubo method All the systems have around 800 atoms. The NEMD methods run for short times; the G-K and Einstein systems need to run longer to generate good statistics. The scripts were all run on a single processor. They all run in a minute or so and produce the accompanying log files and profile files (for velocity or momemtum flux). See the Movies page of the LAMMPS web site (http://lammps.sandia.gov/movies.html), for animations of the NEMD scripts, created using the dump image command. The state point of the LJ fluid is rho* = 0.6, T* = 1.0, and Rcut = 2.5 sigma. This system should have a shear viscosity of about 1.0. ------------- Here is how to extract viscosity from the log file output for each method. The NEMD methods use the formula eta = - dM / Velocity-gradient, where dM = momentum flux in the y-direction, and Vel gradient = dVelX / dY = the change in x-velocity over a distance dY in the y-direction. (1) in.wall.2d mom flux = pxy dVelX = Srate = 2.7 dY = Y box length = 41.99 eta = 0.946 = running average output as last log file column (2) in.nemd.2d mom flux = pxy dVelX = velocity of top box edge = Srate = 2.7 dY = Y box length = 36.51 eta = 1.18 = running average output as last log file column (3) in.mp.2d mom flux = dMom in Y / 2 / Area-perp-to-Y / dTime dMom = -1370.2 from log file, tallied by MP factor of 2 since system is periodic and dMom goes 2 ways Area for 2d = lx dTime = elapsed time in tau for accumulating dMom dVelX = 4th column of log output, from fix ave/spatial dY = 1/2 of Y box length eta = 0.997 = running average output as last log file column (4) in.gk.2d eta is computed directly within the script, by performing a time integration of the formula discussed in Section 6.21 of the manual, -analagous to the formula for thermal conductivity given on the compute +analogous to the formula for thermal conductivity given on the compute heat/flux doc page - the resulting value prints at the end of the run and is in the log file eta = 1.07 (5) in.einstein.2d eta is computed directly within the script, by performing a time integration of the formula discussed in Section 6.21 of the manual, -analagous to the formula for thermal conductivity given on the compute +analogous to the formula for thermal conductivity given on the compute heat/flux doc page - the resulting value prints at the end of the run and is in the log file eta = 1.07