restore gaussian flow example that was lost. tweak input to make it usable for comparing

The input script in.GD is an example simulation using Gaussian dynamics (GD).

The simulation is of a simple 2d Lennard-Jones fluid flowing through a pipe.

For details see online LAMMPS documentation and

Strong and Eaves, J. Phys. Chem. Lett. 7(10) 2016, p. 1907.

Note that the run times and box size are chosen to allow a fast example run.

They are not adequate for a real simulation.

The script has the following parts:

1) initialize variables

    These can be modified to customize the simulation. Note that if the

    pipe dimensions L or d are changed, the geometry should be checked

    by visualizing the coordinates in all.init.lammpstrj.

2) create box

3) set up potential

4) create atoms

5) set up profile-unbiased thermostat (PUT)

    see Evans and Morriss, Phys. Rev. Lett. 56(20) 1986, p. 2172

    By default, this uses boxes which contain on average 8 molecules.

6) equilibrate without GD

7) initialize the center-of-mass velocity and run to achieve steady-state

    The system is initialized with a uniform velocity profile, which

    relaxes over the course of the simulation.

8) collect data

    The data is output in several files:

    GD.out contains the force that GD applies, and the flux in the x- and

        y- directions. The output Jx should be equal to the value of

        J set in section 1, which is 0.1 by default.

    x_profiles contains the velocity, density, and pressure profiles in

        the x-direction. The pressure profile is given by

        (-1/2V)*(c_spa[1] + c_spa[2]), where V is the volume of a

        slice. The pressure profile is computed with IK1, see

        Todd, Evans, and Davis, Phys. Rev. E 52(2) 1995, p. 1627.

        Note that to compare with the pump method, or to

        compute a pressure drop, you must correct this pressure

        profile as described in Strong 2016 above.

    Vy_profile is the velocity profile inside the pipe along the

        y-direction, u_x(y).

#LAMMPS input script

#in.GD

#see README for details

###############################################################################

#initialize variables

clear

#frequency for outputting info (timesteps)

variable        dump_rate       equal 50

variable        thermo_rate     equal 10

#equilibration time (timesteps)

variable        equil           equal 1000

#stabilization time (timesteps to reach steady-state)

variable        stabil          equal 1000

#data collection time (timesteps)

variable        run             equal 2000

#length of pipe

variable        L               equal 30

#width of pipe

variable        d               equal 20

#flux (mass/sigma*tau)

variable        J               equal 0.1

#simulation box dimensions

variable        Lx              equal 100

variable        Ly              equal 40

#bulk fluid density

variable        dens            equal 0.8

#lattice spacing for wall atoms

variable        aWall           equal 1.0 #1.7472

#timestep

variable        ts              equal 0.001

#temperature

variable        T               equal 2.0

#thermostat damping constant

variable        tdamp           equal ${ts}*100

units           lj

dimension       2

atom_style      atomic

###############################################################################

#create box

#create lattice with the spacing aWall

variable        rhoWall         equal ${aWall}^(-2)

lattice         sq ${rhoWall}

#modify input dimensions to be multiples of aWall

variable        L1 equal round($L/${aWall})*${aWall}

variable        d1 equal round($d/${aWall})*${aWall}

variable        Ly1 equal round(${Ly}/${aWall})*${aWall}

variable        Lx1 equal round(${Lx}/${aWall})*${aWall}

#create simulation box

variable        lx2 equal ${Lx1}/2

variable        ly2 equal ${Ly1}/2

region          simbox block -${lx2} ${lx2} -${ly2} ${ly2} 0 0.1 units box

create_box      2 simbox

#####################################################################

#set up potential

mass            1 1.0           #fluid atoms

mass            2 1.0           #wall atoms

pair_style      lj/cut 2.5

pair_modify     shift yes

pair_coeff      1 1 1.0 1.0 2.5

pair_coeff      1 2 1.0 1.0 1.12246

pair_coeff      2 2 0.0 0.0

neigh_modify  exclude type 2 2

timestep        ${ts}

#####################################################################

#create atoms

#create wall atoms everywhere

create_atoms    2 box

#define region which is "walled off"

variable        dhalf equal ${d1}/2

variable        Lhalf equal ${L1}/2

region          walltop block -${Lhalf} ${Lhalf} ${dhalf} EDGE -0.1 0.1 &

                units box

region          wallbot block -${Lhalf} ${Lhalf} EDGE -${dhalf} -0.1 0.1 &

                units box

region          outsidewall union 2 walltop wallbot side out

#remove wall atoms outside wall region

group           outside region outsidewall

delete_atoms    group outside

#remove wall atoms that aren't on edge of wall region

variable        x1 equal ${Lhalf}-${aWall}

variable        y1 equal ${dhalf}+${aWall}

region          insideTop block -${x1} ${x1} ${y1} EDGE -0.1 0.1 units box

region          insideBot block -${x1} ${x1} EDGE -${y1} -0.1 0.1 units box

region          insideWall union 2 insideTop insideBot

group           insideWall region insideWall

delete_atoms    group insideWall

#define new lattice, to give correct fluid density

#y lattice const must be a multiple of aWall

variable        atrue equal ${dens}^(-1/2)

variable        ay equal round(${atrue}/${aWall})*${aWall}

#choose x lattice const to give correct density

variable        ax equal (${ay}*${dens})^(-1)

#change Lx to be multiple of ax

variable        Lx1 equal round(${Lx}/${ax})*${ax}

variable        lx2 equal ${Lx1}/2

change_box      all x final -${lx2} ${lx2} units box

#define new lattice

lattice         custom ${dens} &

                a1 ${ax} 0.0 0.0 a2 0.0 ${ay} 0.0 a3 0.0 0.0 1.0 &

                basis 0.0 0.0 0.0

#fill in rest of box with bulk particles

variable        delta equal 0.001

variable        Ldelt equal ${Lhalf}+${delta}

variable        dDelt equal ${dhalf}-${delta}

region          left block EDGE -${Ldelt} EDGE EDGE -0.1 0.1 units box

region          right block ${Ldelt} EDGE EDGE EDGE -0.1 0.1 units box

region          pipe block -${Ldelt} ${Ldelt} -${dDelt} ${dDelt} -0.1 0.1 &

                units box

region          bulk union 3 left pipe right

create_atoms    1 region bulk

group           bulk type 1

group           wall type 2

#remove atoms that are too close to wall

delete_atoms    overlap 0.9 bulk wall

neighbor        0.3 bin

neigh_modify    delay 0 every 1 check yes

neigh_modify    exclude group wall wall

velocity        bulk create $T 78915 dist gaussian rot yes mom yes loop geom

#####################################################################

#set up PUT

#see Evans and Morriss, Phys. Rev. Lett. 56(20) 1986, p. 2172

#average number of particles per box, Evans and Morriss used 2.0

variable        NperBox equal 8.0

#calculate box sizes

variable        boxSide equal sqrt(${NperBox}/${dens})

variable        nX equal round(lx/${boxSide})

variable        nY equal round(ly/${boxSide})

variable        dX equal lx/${nX}

variable        dY equal ly/${nY}

#temperature of fluid (excluding wall)

compute         myT bulk temp

#profile-unbiased temperature of fluid

compute         myTp bulk temp/profile 1 1 0 xy ${nX} ${nY}

#thermo setup

thermo          ${thermo_rate}

thermo_style    custom step c_myT c_myTp etotal press

#dump initial configuration

# dump            55 all custom 1 all.init.lammpstrj id type x y z vx vy vz

# dump            56 wall custom 1 wall.init.lammpstrj id type x y z

# dump_modify     55 sort id

# dump_modify     56 sort id

run             0

# undump          55

# undump          56

#####################################################################

#equilibrate without GD

fix             nvt bulk nvt temp $T $T ${tdamp}

fix_modify      nvt temp myTp

fix             2 bulk enforce2d

run             ${equil}

#####################################################################

#initialize the COM velocity and run to achieve steady-state

#calculate velocity to add: V=J/rho_total

variable        Vadd            equal $J*lx*ly/count(bulk)

#first remove any COM velocity, then add back the streaming velocity

velocity        bulk zero linear

velocity        bulk set ${Vadd} 0.0 0.0 units box sum yes mom no

fix             GD bulk flow/gauss 1 0 0 #energy yes

#fix_modify     GD energy yes

run             ${stabil}

#####################################################################

#collect data

#print the applied force and total flux to ensure conservation of Jx

variable        Fapp equal f_GD[1]

compute         vxBulk bulk reduce sum vx

compute         vyBulk bulk reduce sum vy

variable        invVol equal 1.0/(lx*ly)

variable        jx equal c_vxBulk*${invVol}

variable        jy equal c_vyBulk*${invVol}

variable        curr_step equal step

variable        p_Fapp format Fapp %.3f

variable        p_jx format jx %.5g

variable        p_jy format jy %.5g

fix             print_vCOM all print ${dump_rate} &

                "${curr_step} ${p_Fapp} ${p_jx} ${p_jy}" file GD.out screen no &

                title "timestep Fapp Jx Jy"

#compute IK1 pressure profile

#see Todd, Evans, and Davis, Phys. Rev. E 52(2) 1995, p. 1627

#use profile-unbiased temperature to remove the streaming velocity

#from the kinetic part of the pressure

compute         spa bulk stress/atom myTp

#for the pressure profile, use the same grid as the PUT

compute         chunkX bulk chunk/atom bin/1d x lower ${dX} units box

#output pressure profile and other profiles

#the pressure profile is (-1/2V)*(c_spa[1] + c_spa[2]), where

#V is the volume of a slice

fix             profiles bulk ave/chunk 1 1 ${dump_rate} chunkX &

                vx density/mass  c_spa[1] c_spa[2] &

                file x_profiles ave running overwrite

#compute velocity profile across the pipe with a finer grid

variable        dYnew equal ${dY}/10

compute         chunkY bulk chunk/atom bin/1d y center ${dYnew} units box &

                region pipe

fix             velYprof bulk ave/chunk 1 1 ${dump_rate} chunkY &

                vx file Vy_profile ave running overwrite

#full trajectory

# dump          7 bulk custom ${dump_rate} bulk.lammpstrj id type x y z

# dump_modify     7 sort id

run             ${run}

timestep Fapp Jx Jy

2050 -215.835 0.1 -0.002562

2100 -220.455 0.1 -0.0019705

2150 55.212 0.1 -0.0028338

2200 87.052 0.1 -0.0042335

2250 -62.998 0.1 -0.0045646

2300 71.630 0.1 -0.0039858

2350 43.159 0.1 -0.0029771

2400 109.930 0.1 -0.0018522

2450 110.735 0.1 -0.0011188

2500 107.071 0.1 0.0005978

2550 335.449 0.1 0.0010164

2600 159.694 0.1 -0.00015953

2650 6.532 0.1 -0.0004907

2700 65.524 0.1 -0.00093116

2750 79.662 0.1 -0.0033425

2800 69.846 0.1 -0.0055377

2850 122.175 0.1 -0.00721

2900 32.456 0.1 -0.0086166

2950 -85.137 0.1 -0.01107

3000 154.735 0.1 -0.011337

3050 72.979 0.1 -0.0095316

3100 -24.457 0.1 -0.0098708

3150 -0.383 0.1 -0.0094961

3200 132.434 0.1 -0.011524

3250 48.222 0.1 -0.014966

3300 -73.186 0.1 -0.016999

3350 172.062 0.1 -0.018554

3400 106.144 0.1 -0.021202

3450 -22.860 0.1 -0.01949

3500 22.120 0.1 -0.016033

3550 -254.920 0.1 -0.012172

3600 -147.218 0.1 -0.011162

3650 -12.508 0.1 -0.010255

3700 81.846 0.1 -0.0085117

3750 -79.406 0.1 -0.0061294

3800 -34.994 0.1 -0.0026239

3850 94.992 0.1 -0.0015312

3900 -0.345 0.1 -0.0011157

3950 -88.693 0.1 -0.0018929

4000 156.029 0.1 -0.0024547