Update energy and time scaling in bash scripts

tut0_lj_fluid/init.sh

@@ -27,9 +27,9 @@ rc_0=1.122462048309373 # 2^(1/6) - the cutoff distance for a purely
 
 sigma=1.0 # Set all atoms to have the same size
 
-# Rescale epsilon such that the minimum of the truncated LJ potential
-# actually reaches -epsilon
-function get_norm() {
+# Rescale epsilon such that the minimum of the truncated and shifted LJ
+# potential actually reaches -epsilon
+function rescale_energy() {
     local e=$1
     local rc=$2
     local s=$3
@@ -37,22 +37,29 @@ function get_norm() {
 norm = 1.0+4.0*(($s/$rc)**12.0-($s/$rc)**6.0)
 print($e/norm if norm > 0.0 else $e)")
 }
-epsilon_11=$(get_norm $epsilon_11 $rc_11 $sigma)
-epsilon_12=$(get_norm $epsilon_12 $rc_12 $sigma)
-epsilon_22=$(get_norm $epsilon_22 $rc_22 $sigma)
-
-# Dump frequency and length of simulation
-dt=0.01 # Set timestep size in LJ simulation time unit
-dump_printfreq=100 # In LJ simulation time unit
-thermo_printfreq=100
-run_init_time=100
-run_time=10000
-
-# Convert from simulation time unit to timesteps
-dump_printfreq=$($pyx -c "print(int($dump_printfreq/$dt))")
-thermo_printfreq=$($pyx -c "print(int($thermo_printfreq/$dt))")
-run_init_time=$($pyx -c "print(int($run_init_time/$dt))")
-run_time=$($pyx -c "print(int($run_time/$dt))")
+epsilon_11=$(rescale_energy $epsilon_11 $rc_11 $sigma)
+epsilon_12=$(rescale_energy $epsilon_12 $rc_12 $sigma)
+epsilon_22=$(rescale_energy $epsilon_22 $rc_22 $sigma)
+
+# Simulation run times (in simulation time unit)
+dt=0.01 # Size of each timestep
+run_init_time=100    # Equilibration run time
+run_time=10000       # Main simulation run time
+
+# Dump frequencies (in simulation time unit)
+dump_printfreq=100   # Frequency for dumping atom positions
+thermo_printfreq=100 # Frequency for dumping thermodynamics info
+
+# Convert run times and dump frequencies to be in timesteps
+function get_timestep() {
+    local tau=$1
+    local dt=$2
+    echo $($pyx -c "print(int($tau/$dt))")
+}
+run_init_time=$(get_timestep $run_init_time $dt)
+run_time=$(get_timestep $run_time $dt)
+dump_printfreq=$(get_timestep $dump_printfreq $dt)
+thermo_printfreq=$(get_timestep $thermo_printfreq $dt)
 
 # Set the box boundaries
 xlo=$($pyx -c "print(-int(${lx}/2.0))")
@@ -66,10 +73,8 @@ zhi=$($pyx -c "print(int(${lz}/2.0))")
 pos_equil_file="pos_equil.lammpstrj"
 pos_file="pos.lammpstrj"
 
-# Required programs and files
-atom_py="create_atoms.py"
-
 # Generate the LAMMPS input file
+atom_py="create_atoms.py"
 traj_in="traj.in"
 $pyx $atom_py $natoms $ntypes $lx $ly $lz $seed $traj_in
 

tut1_polymer/init.sh

@@ -31,9 +31,9 @@ R_0=1.6
 # Angle potential
 l_p=3.0 # Set the persistence length of the polymer
 
-# Rescale epsilon such that the minimum of the truncated LJ potential
-# actually reaches -epsilon
-function get_norm() {
+# Rescale epsilon such that the minimum of the truncated and shifted LJ
+# potential actually reaches -epsilon
+function rescale_energy() {
     local e=$1
     local rc=$2
     local s=$3
@@ -41,22 +41,29 @@ function get_norm() {
 norm = 1.0+4.0*(($s/$rc)**12.0-($s/$rc)**6.0)
 print($e/norm if norm > 0.0 else $e)")
 }
-epsilon=$(get_norm $epsilon $rc $sigma)
-
-# Dump frequency and length of simulation
-dt=0.01 # Set timestep size in LJ simulation time unit
-dump_printfreq=100 # In LJ simulation time unit
-thermo_printfreq=100
-run_init_time_1=100
-run_init_time_2=100
-run_time=10000
-
-# Convert from simulation time unit to timesteps
-dump_printfreq=$($pyx -c "print(int($dump_printfreq/$dt))")
-thermo_printfreq=$($pyx -c "print(int($thermo_printfreq/$dt))")
-run_init_time_1=$($pyx -c "print(int($run_init_time_1/$dt))")
-run_init_time_2=$($pyx -c "print(int($run_init_time_2/$dt))")
-run_time=$($pyx -c "print(int($run_time/$dt))")
+epsilon=$(rescale_energy $epsilon $rc $sigma)
+
+# Simulation run times (in simulation time unit)
+dt=0.01 # Size of each timestep
+run_init_time_1=100  # Equilibration time with harmonic bonds
+run_init_time_2=100  # Equilibration time with FENE bonds
+run_time=10000       # Main simulation run time
+
+# Dump frequencies (in simulation time unit)
+dump_printfreq=100   # Frequency for dumping bead positions
+thermo_printfreq=100 # Frequency for dumping thermodynamics info
+
+# Convert run times and dump frequencies to be in timesteps
+function get_timestep() {
+    local tau=$1
+    local dt=$2
+    echo $($pyx -c "print(int($tau/$dt))")
+}
+run_init_time_1=$(get_timestep $run_init_time_1 $dt)
+run_init_time_2=$(get_timestep $run_init_time_2 $dt)
+run_time=$(get_timestep $run_time $dt)
+dump_printfreq=$(get_timestep $dump_printfreq $dt)
+thermo_printfreq=$(get_timestep $thermo_printfreq $dt)
 
 # Set the box boundaries
 xlo=$($pyx -c "print(-int(${lx}/2.0))")
@@ -70,10 +77,8 @@ zhi=$($pyx -c "print(int(${lz}/2.0))")
 pos_equil_file="pos_equil.lammpstrj"
 pos_file="pos.lammpstrj"
 
-# Required programs and files
-polymer_py="create_polymer.py"
-
 # Generate the LAMMPS input file
+polymer_py="create_polymer.py"
 traj_in="traj.in"
 $pyx $polymer_py $nbeads $ntypes $sigma $lx $ly $lz $seed $traj_in
 
@@ -150,7 +155,7 @@ read_data $(basename $init_file)
 # screen and/or log file
 thermo ${thermo_printfreq} # Set the output frequency
 compute gyr all gyration # Compute the radius of gyration of the polymer
-thermo_style custom step temp epair emol vol c_gyr
+thermo_style custom step temp epair vol c_gyr
 
 # Set up a dump to output bead positions
 dump 1 all custom ${dump_printfreq} $(basename $pos_equil_file) &

tut2_dna_protein/init.sh

@@ -31,9 +31,9 @@ R_0=1.6
 # Angle potential
 l_p=3.0 # Set the persistence length of the DNA/chromatin fibre
 
-# Rescale epsilon such that the minimum of the truncated LJ potential
-# actually reaches -epsilon
-function get_norm() {
+# Rescale epsilon such that the minimum of the truncated and shifted LJ
+# potential actually reaches -epsilon
+function rescale_energy() {
     local e=$1
     local rc=$2
     local s=$3
@@ -41,22 +41,29 @@ function get_norm() {
 norm = 1.0+4.0*(($s/$rc)**12.0-($s/$rc)**6.0)
 print($e/norm if norm > 0.0 else $e)")
 }
-epsilon=$(get_norm $epsilon $rc $sigma)
-
-# Dump frequency and length of simulation
-dt=0.01 # Set timestep size in LJ simulation time unit
-dump_printfreq=100 # In LJ simulation time unit
-thermo_printfreq=100
-run_init_time_1=100
-run_init_time_2=100
-run_time=10000
-
-# Convert from simulation time unit to timesteps
-dump_printfreq=$($pyx -c "print(int($dump_printfreq/$dt))")
-thermo_printfreq=$($pyx -c "print(int($thermo_printfreq/$dt))")
-run_init_time_1=$($pyx -c "print(int($run_init_time_1/$dt))")
-run_init_time_2=$($pyx -c "print(int($run_init_time_2/$dt))")
-run_time=$($pyx -c "print(int($run_time/$dt))")
+epsilon=$(rescale_energy $epsilon $rc $sigma)
+
+# Simulation run times (in simulation time unit)
+dt=0.01 # Size of each timestep
+run_init_time_1=100  # Equilibration time with harmonic bonds
+run_init_time_2=100  # Equilibration time with FENE bonds
+run_time=10000       # Main simulation run time
+
+# Dump frequencies (in simulation time unit)
+dump_printfreq=100   # Frequency for dumping bead positions
+thermo_printfreq=100 # Frequency for dumping thermodynamics info
+
+# Convert run times and dump frequencies to be in timesteps
+function get_timestep() {
+    local tau=$1
+    local dt=$2
+    echo $($pyx -c "print(int($tau/$dt))")
+}
+run_init_time_1=$(get_timestep $run_init_time_1 $dt)
+run_init_time_2=$(get_timestep $run_init_time_2 $dt)
+run_time=$(get_timestep $run_time $dt)
+dump_printfreq=$(get_timestep $dump_printfreq $dt)
+thermo_printfreq=$(get_timestep $thermo_printfreq $dt)
 
 # Set the box boundaries
 xlo=$($pyx -c "print(-int(${lx}/2.0))")
@@ -70,10 +77,8 @@ zhi=$($pyx -c "print(int(${lz}/2.0))")
 pos_equil_file="pos_equil.lammpstrj"
 pos_file="pos.lammpstrj"
 
-# Required programs and files
-dna_py="create_dna_protein.py"
-
 # Generate the LAMMPS input file
+dna_py="create_dna_protein.py"
 traj_in="traj.in"
 $pyx $dna_py $nbeads $nprots $sigma $lx $ly $lz $seed $traj_in
 
@@ -155,7 +160,7 @@ group poly type 1
 # Set the params for dumping thermodynamic properties
 thermo ${thermo_printfreq}
 compute gyr poly gyration
-thermo_style custom step temp epair emol vol c_gyr
+thermo_style custom step temp epair vol c_gyr
 
 # Set up a dump to output bead positions
 dump 1 all custom ${dump_printfreq} $(basename $pos_equil_file) &