updated fix_langevin and made example

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

#Atoms types - mass - charge#

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

#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#

variable Ar equal 1

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

#Atom Masses#

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

mass ${Ar}      39.903

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

#Pair Potentials - Tersoff#

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

pair_style      lj/cubic

pair_coeff      * * 0.0102701 3.42

###############################mm

# Atom style - charge/vdw/bonded#

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

atom_style full

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

#Units Metal : eV       - ps - angstrom - bar#

#      Real  : kcal/mol - fs - angstrom - atm#

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

units   	metal

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

#Run number#

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

variable run_no   equal   0             # is it a restart?

variable res_no   equal   ${run_no}-1   # restart file number 

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

#Random Seeds and Domain Decomposition#

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

variable iseed0 equal 2357  

variable iseed1 equal 26488 

variable iseed2 equal 10669 

processors      * * *

###########

#Data File#

###########

variable inpfile   string   argon.lmp

variable resfile   string   final_restart.${res_no}

variable ff_file   string   ff-argon.lmp

##########

#Run Type#

##########

variable minimise equal    0   #Energy Minimization

variable md       equal    1   #Plain MD

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

#Molecular Dynamics Parameters#

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

variable run_no   equal        0       # is it a restart?

variable ens      equal        9       # ensemble (0=nve, 1=nvt, 2=npt, 3=ber, 4=lang, 5=stoc, 6=vres, 7=stoch, 8=gjf)

variable ts       equal        0.120   # simulation timestep (time units)

variable nequil   equal        0       # number of equilibration steps

variable nsteps   equal   200000       # number of MD steps

#variable nsteps   equal   20       # number of MD steps

variable temp_s   equal       10       # starting temperature

variable temp_f   equal       10       # final simulation temperature

variable trel     equal        1       # thermostat relaxation time

variable tscale   equal        1       # thermostat relaxation freq - vel rescaling only

variable deltat   equal        1       # maximum temperature change - vel rescaling only

variable npttype  string       iso     # type of NPT (iso, aniso, tri, z...)

variable pres     equal        1.01325 # pressure (NPT runs only)

variable prel     equal        1.0     # barostat relaxation time

neighbor 1 bin

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

#Output Parameters#

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

variable ntraj    equal     1000       # trajectory output frequency - all system

variable ntraj_s  equal     -100       # trajectory output frequency - solute only

variable nthermo  equal      200       # thermodynamic data output frequency

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

#Energy Minimization Parameters#

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

variable mtraj    equal        1       # trajectory output frequency - all system

variable etol     equal     1e-5       # % change in energy

variable ftol     equal     1e-5       # max force threshold (force units)

variable maxiter  equal    10000       # max # of iterations

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

#3D Periodic Simulation#

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

boundary 	p p p

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

#Reading the input structure#

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

if "${run_no} == 0" then "read_data ${inpfile}" else "read_restart ${resfile}"

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

#Force Field#

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

include ${ff_file}

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

#Thermodynamic Output#

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

variable str_basic string 'step time pe temp press'

#MD ensemble (0=nve, 1=nvt, 2=npt, 3=ber, 4=lang, 5=stoc, 6=vres)

variable str_ens string ' '

if "${ens} == 0" then "variable str_ens string 'etotal'"

if "${ens} == 2" then "variable str_ens string 'vol pxx pyy pzz cella cellb cellc cellakpha cellbeta cellgamma'"

#Variable for a gulp friend output

if "${ens} >= 0" then "thermo_style custom time temp  pe  etotal press vol  cpu" &

        "thermo ${nthermo}" &

        "thermo_modify flush yes" 

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

#Energy Minimization#

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

if "${minimise} <= 0 || ${run_no} > 0" then "jump SELF end_minimise"

  print "Doing CG minimisation"

  dump mdcd all dcd ${mtraj} min.dcd

  dump_modify mdcd unwrap yes

  min_style cg

  min_modify line quadratic

  minimize ${etol} ${ftol} ${maxiter} ${maxiter}

  reset_timestep 0

  undump mdcd

label end_minimise

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

#Timestep in ps#

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

timestep ${ts}

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

#Restart file#

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

restart 100000 restart.1 restart.2

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

#Trajectory output#

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

#dump  xyz all atom 1000 silicon.lammpstrj

if "${ntraj} > 0" then &

  "dump 1 all dcd ${ntraj} trajectory.${run_no}.dcd" &

  "dump_modify 1 unwrap yes"

fix mom all momentum 1 linear 1 1 1

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

#Ensembles (0=nve,1=nvt, 2=npt, 3=ber, 4=lang, 5=stoc, 6=vres)#

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

if "${md} > 0" then 'print "Setting up the ensembles"' &

  'if "${run_no} == 0" then "velocity all create ${temp_s} ${iseed0} mom yes dist gaussian"' &

  'if "${ens} == 0" then "fix nve all nve"' &

  'if "${ens} == 1" then "fix nvt all nvt temp ${temp_s} ${temp_f} ${trel} tchain 5"' &

  'if "${ens} == 2" then "fix npt all npt temp ${temp_s} ${temp_f} ${trel} ${npttype} ${pres} ${pres} ${prel} tchain 5 pchain 5 mtk yes"' &

  'if "${ens} == 3" then "fix nve all nve" "fix ber all temp/berendsen ${temp_s} ${temp_f} ${trel}"' &

  'if "${ens} == 4" then "fix nve all nve" "fix lang all langevin ${temp_s} ${temp_f} ${trel} ${iseed1} tally yes zero yes"' &

  'if "${ens} == 5" then "fix nve all nve" "fix stoch all temp/csvr ${temp_s} ${temp_f} ${trel} ${iseed1}"' &

  'if "${ens} == 6" then "fix nve all nve" "fix stoch all temp/csld ${temp_s} ${temp_f} ${trel} ${iseed1}"' &

  'if "${ens} == 7" then "fix nve all nve" "fix vres all temp/rescale ${tscale} ${temp_s} ${temp_f} ${tmin} ${tmax}"' &

  'if "${ens} == 8" then "fix nve all nve" "fix lang all langevin ${temp_s} ${temp_f} ${trel} ${iseed1} gjf yes"' &

  'if "${ens} == 9" then "fix nve all nve" "fix lang all langevin ${temp_s} ${temp_f} ${trel} ${iseed1} gjf yes halfstep yes"'

if "${md} > 0" then "print 'Doing Molecular dynamics'" &

        "run ${nsteps}" &

        "write_restart final_restart.${run_no}"