# Rodrigo Freitas - rodrigohb@gmail.com

#-------------------------- Initial Setup -------------------------------------#

  units             metal

  atom_modify       map array 

#------------------------------------------------------------------------------#

#-------------------------- Atomic Setup --------------------------------------#

  lattice           sc 10.0

  region            my_region block 0 8 0 8 0 8

  create_box        1 my_region

  create_atoms      1 box

  mass              * 50.0

  atom_modify       sort 0 0

#------------------------------------------------------------------------------#

#---------------------- Thermostat & Barostat ---------------------------------#

  fix               f1 all nve

  fix               f2 all spring/self 10.0

  fix               f3 all ti/rs 1.0 0.5 2000 1000 function 2

  fix               f4 all langevin 100.0 100.0 0.5 1230917 zero yes

#------------------------------------------------------------------------------#

#------------------ Computes, variables & modifications -----------------------#

  compute           Tcm all temp/com

  fix_modify        f4 temp Tcm

  variable          step    equal step

  variable          U       equal f_f2/atoms

  variable          lambda  equal f_f3[1]

  variable          dlambda equal f_f3[2]

#------------------------------------------------------------------------------#

#------------------------- Thermo stuff ---------------------------------------#

  thermo            0

  timestep          0.001

#------------------------------------------------------------------------------#

#------------------------- Running the Simulation -----------------------------#

  velocity          all create 200.0 8128481 mom yes rot yes dist gaussian

  # Forward. 

  run               999

  fix               f5 all print 1 "${step} ${U} ${lambda} ${dlambda}" & 

                    screen no file forward_rs.dat &

                    title '# step U lambda dlambda'

  run               2001

  unfix             f5

  # Backward.

  run               999

  fix               f5 all print 1 "${step} ${U} ${lambda} ${dlambda}" & 

                    screen no file backward_rs.dat &

                    title '# step U lambda dlmabda'

  run               2001

  unfix             f5

#------------------------------------------------------------------------------#

# Rodrigo Freitas - rodrigohb@gmail.com

# Rodrigo Freitas - rodrigof@berkeley.edu

#

# Description: nonequilibrium thermodynamic integration. Further details in:

# R. Freitas, M. Asta, and M. de Koning, Computational Materials Science, (2016)

# http://dx.doi.org/10.1016/j.commatsci.2015.10.050

#-------------------------- Initial Setup -------------------------------------#

#--------------------------- System setup -------------------------------------#

  units             metal

  atom_modify       map array 

#------------------------------------------------------------------------------#

#-------------------------- Atomic Setup --------------------------------------#

  lattice           sc 10.0

  region            my_region block 0 8 0 8 0 8

  create_box        1 my_region

  lattice           fcc 3.621

  region            sim_box block 0 4 0 4 0 4

  create_box        1 sim_box

  create_atoms      1 box

  mass              * 50.0

  atom_modify       sort 0 0

  pair_style        eam/alloy

  pair_coeff        * * ../../../../potentials/Cu_mishin1.eam.alloy Cu

#------------------------------------------------------------------------------#

#---------------------- Thermostat & Barostat ---------------------------------#

#---------------------- Simulation setup --------------------------------------#

  # The order of the fix commands is important: thermostat AFTER thermodynamic integration force interpolation.

  fix               f1 all nve

  fix               f2 all spring/self 10.0

  fix               f3 all ti/spring 50.0 2000 1000 function 2

  fix               f4 all langevin 100.0 100.0 0.5 1230917 zero yes

#------------------------------------------------------------------------------#

  fix               f2 all ti/spring 3.728 2000 1000 function 2

  fix               f3 all langevin 100.0 100.0 0.1 666 zero yes

  # Compute temperature using center-of-mass coordinates.

  compute           c1 all temp/com

  fix_modify        f3 temp c1

#------------------ Computes, variables & modifications -----------------------#

  compute           Tcm all temp/com

  fix_modify        f4 temp Tcm

  # Output variables.

  variable          step    equal step

  variable          dU      equal (f_f2-f_f3)/atoms

  variable          lambda  equal f_f3[1]

  variable          dlambda equal f_f3[2]

#------------------------------------------------------------------------------#

  variable          dU      equal pe-f_f2 # Used on the integration.

  variable          lambda  equal f_f2[1] # For reference or integration.

  variable          dlambda equal f_f2[2] # For reference or integration.

#------------------------- Thermo stuff ---------------------------------------#

  thermo            0

  timestep          0.001

  # Thermo output.

  thermo_style      custom step pe

  thermo            100

#------------------------------------------------------------------------------#

#------------------------- Running the Simulation -----------------------------#

  velocity          all create 200.0 8128481 mom yes rot yes dist gaussian

#------------------------- Running the simulation -----------------------------#

  # Setup initial velocities to accelerate equilibration.

  velocity          all create 200.0 9999 mom yes rot yes dist gaussian

  # Forward. 

  run               999

  fix               f5 all print 1 "${step} ${dU} ${lambda} ${dlambda}" & 

                    screen no file forward_fl.dat &

                    title '# step dU lambda dlambda'

  run               2001

  unfix             f5

  # Forward nonequilibrium thermodynamic integration.

  run               1000 # Equilibrate system at lambda = 0.

  fix               f4 all print 1 "${step} ${dU} ${lambda} ${dlambda}" & 

                    screen no file forward.dat title '# step dU lambda dlambda'

  run               2000 # Lambda switching: 0 -> 1.

  unfix             f4

  # Backward.

  run               999

  fix               f5 all print 1 "${step} ${dU} ${lambda} ${dlambda}" & 

                    screen no file backward_fl.dat &

                    title '# step dU lambda dlambda'

  run               2001

  unfix             f5

  # Backward nonequilibrium thermodynamic integration.

  run               1000 # Equilibrate system at lambda = 1.

  fix               f4 all print 1 "${step} ${dU} ${lambda} ${dlambda}" & 

                    screen no file backward.dat title '# step dU lambda dlambda'

  run               2000 # Lambda switching: 1 -> 0.

  unfix             f4

#------------------------------------------------------------------------------#