diff --git a/calphy/alchemy.py b/calphy/alchemy.py
index 27d1be1..2ac405c 100644
--- a/calphy/alchemy.py
+++ b/calphy/alchemy.py
@@ -112,9 +112,8 @@ def run_averaging(self):
         self.check_if_melted(lmp, "traj.equilibration_stage2.dat")
 
         #close object and process traj
-        lmp = ph.write_data(lmp, "current.data")
+        lmp = ph.write_data(lmp, "conf.equilibration.data")
         lmp.close()
-        self.process_traj("traj.equilibration_stage2.dat", "conf.equilibration.data")
 
 
     
diff --git a/calphy/clitools.py b/calphy/clitools.py
new file mode 100644
index 0000000..199bea6
--- /dev/null
+++ b/calphy/clitools.py
@@ -0,0 +1,92 @@
+import os
+import numpy as np
+import shutil
+import argparse as ap
+import subprocess
+import yaml
+import time
+import datetime
+
+from calphy.input import read_inputfile, load_job, save_job
+from calphy.liquid import Liquid
+from calphy.solid import Solid
+from calphy.alchemy import Alchemy
+
+def _generate_job(calc, simfolder):
+    if calc.mode == "alchemy" or calc.mode == "composition_scaling":
+        job = Alchemy(calculation=calc, simfolder=simfolder)
+        return job
+    else:
+        if calc.reference_phase == "liquid":
+            job = Liquid(calculation=calc, simfolder=simfolder)
+            return job
+        else:
+            job = Solid(calculation=calc, simfolder=simfolder)
+            return job
+
+
+def run_averaging():
+    arg = ap.ArgumentParser()
+    arg.add_argument("-i", "--input", required=True, type=str,
+    help="name of the input file")
+    arg.add_argument("-k", "--kernel", required=True, type=int, 
+    help="kernel number of the calculation to be run.")
+    args = vars(arg.parse_args())
+    kernel = args["kernel"]
+    calculations = read_inputfile(args["input"])
+    calc = calculations[kernel]
+
+    simfolder = calc.create_folders()
+    job = _generate_job(calc, simfolder)
+    os.chdir(simfolder)
+    
+
+    job.run_averaging()
+    save_job(job)
+
+
+def process_averaging():
+    arg = ap.ArgumentParser()
+    arg.add_argument("-i", "--input", required=True, type=str,
+    help="name of the input file")
+    arg.add_argument("-k", "--kernel", required=True, type=int, 
+    help="kernel number of the calculation to be run.")
+    args = vars(arg.parse_args())
+    kernel = args["kernel"]
+    calculations = read_inputfile(args["input"])
+    calc = calculations[kernel]
+
+    job = load_job(calc.savefile)
+    job.process_averaging_results()
+    save_job(job)
+
+def run_integration():
+    arg = ap.ArgumentParser()
+    arg.add_argument("-i", "--input", required=True, type=str,
+    help="name of the input file")
+    arg.add_argument("-k", "--kernel", required=True, type=int, 
+    help="kernel number of the calculation to be run.")
+    args = vars(arg.parse_args())
+    kernel = args["kernel"]
+    calculations = read_inputfile(args["input"])
+    calc = calculations[kernel]
+
+    job = load_job(calc.savefile)
+    job.run_integration()
+    save_job(job)
+
+def process_integration():
+    arg = ap.ArgumentParser()
+    arg.add_argument("-i", "--input", required=True, type=str,
+    help="name of the input file")
+    arg.add_argument("-k", "--kernel", required=True, type=int, 
+    help="kernel number of the calculation to be run.")
+    args = vars(arg.parse_args())
+    kernel = args["kernel"]
+    calculations = read_inputfile(args["input"])
+    calc = calculations[kernel]
+
+    job = load_job(calc.savefile)
+    job.thermodynamic_integration()
+    job.submit_report()
+    save_job(job)
diff --git a/calphy/helpers.py b/calphy/helpers.py
index ca23f15..8a2e6be 100644
--- a/calphy/helpers.py
+++ b/calphy/helpers.py
@@ -32,9 +32,20 @@
 from ase.io import read, write
 from pyscal.trajectory import Trajectory
 
+class LammpsScript:
+    def __init__(self):
+        self.script = []
+
+    def command(self, command_str):
+        self.script.append(command_str)
+
+    def write(self, infile):
+        with open(infile, 'w') as fout:
+            for line in self.script:
+                fout.write(f'{line}\n')
 
 def create_object(cores, directory, timestep, cmdargs=None, 
-    init_commands=None):
+    init_commands=None, script_mode=False):
     """
     Create LAMMPS object
 
@@ -53,9 +64,12 @@ def create_object(cores, directory, timestep, cmdargs=None,
     -------
     lmp : LammpsLibrary object
     """
-    lmp = LammpsLibrary(
-        cores=cores, working_directory=directory, cmdargs=cmdargs
-    )
+    if script_mode:
+        lmp = LammpsScript()
+    else:
+        lmp = LammpsLibrary(
+            cores=cores, working_directory=directory, cmdargs=cmdargs
+        )
 
     commands = [["units", "metal"], 
     ["boundary", "p p p"],
@@ -150,8 +164,10 @@ def set_potential(lmp, options, ghost_elements=0):
     -------
     lmp : LammpsLibrary object
     """
-    lmp.pair_style(options.pair_style_with_options[0])
-    lmp.pair_coeff(options.pair_coeff[0])
+    #lmp.pair_style(options.pair_style_with_options[0])
+    #lmp.pair_coeff(options.pair_coeff[0])
+    lmp.command(f'pair_style {options.pair_style_with_options[0]}')
+    lmp.command(f'pair_coeff {options.pair_coeff[0]}')
 
     lmp = set_mass(lmp, options, ghost_elements=ghost_elements)
 
diff --git a/calphy/input.py b/calphy/input.py
index d46a577..c84af5d 100644
--- a/calphy/input.py
+++ b/calphy/input.py
@@ -166,6 +166,9 @@ def __init__(self):
         self._reference_phase = None
         self._lattice_constant = 0
         self._repeat = [1, 1, 1]
+        self._script_mode = False
+        self._lammps_executable = None
+        self._mpi_executable = None
         self._npt = True
         self._n_equilibration_steps = 25000
         self._n_switching_steps = 50000
@@ -242,20 +245,13 @@ def to_string(self, val):
         else:
             return str(val)
 
-    def to_bool(self, val):        
-        if val == "True":
-            val = True
-        elif val == "False":
-            val = False
-        elif val == 1:
-            val = True
-        elif val == 0:
-            val = False
-        elif val == "1":
-            val = True
-        elif val == "0":
-            val = False
-        return val
+    def to_bool(self, val):
+        if val in ["True", "true", 1, "1", True]:
+            return True
+        elif val in ["False", "false", 0, "0", False]:
+            return False
+        else:
+            raise ValueError(f'Unknown bool input of type {type(val)} with value {val}')
 
     @property
     def element(self):
@@ -531,7 +527,39 @@ def npt(self):
     
     @npt.setter
     def npt(self, val):
-        self._npt = val
+        val = self.to_bool(val)
+        if isinstance(val, bool):
+            self._npt = val
+        else:
+            raise TypeError("NPT should be either True/False")
+
+    @property
+    def script_mode(self):
+        return self._script_mode
+    
+    @script_mode.setter
+    def script_mode(self, val):
+        val = self.to_bool(val)
+        if isinstance(val, bool):
+            self._script_mode = val
+        else:
+            raise TypeError("script mode should be either True/False")
+
+    @property
+    def lammps_executable(self):
+        return self._lammps_executable
+    
+    @lammps_executable.setter
+    def lammps_executable(self, val):
+        self._lammps_executable = val
+
+    @property
+    def mpi_executable(self):
+        return self._mpi_executable
+    
+    @mpi_executable.setter
+    def mpi_executable(self, val):
+        self._mpi_executable = val
     
     @property
     def n_equilibration_steps(self):
@@ -664,7 +692,12 @@ def create_identifier(self):
             identistring = "-".join([self.folder_prefix, prefix, l, str(ts), str(ps)])
         return identistring
 
-    def create_folders(self, prefix=None):
+    def get_folder_name(self):
+        identistring = self.create_identifier()
+        simfolder = os.path.join(os.getcwd(), identistring)
+        return simfolder
+
+    def create_folders(self):
         """
         Create the necessary folder for calculation
 
@@ -678,11 +711,7 @@ def create_folders(self, prefix=None):
         folder : string
             create folder
         """
-        identistring = self.create_identifier()
-        if prefix is None:
-            simfolder = os.path.join(os.getcwd(), identistring)
-        else:
-            simfolder = os.path.join(prefix, identistring)
+        simfolder = self.get_folder_name()
 
         #if folder exists, delete it -> then create
         try:
@@ -706,6 +735,10 @@ def generate(cls, indata):
             calc = cls()
             calc.element = indata["element"]
             calc.mass = indata["mass"]
+            calc.script_mode = indata["script_mode"]
+            calc.lammps_executable = indata["lammps_executable"]
+            calc.mpi_executable = indata["mpi_executable"]
+
             if "md" in indata.keys():
                 calc.md.add_from_dict(indata["md"])
             if "queue" in indata.keys():
@@ -735,6 +768,12 @@ def read_file(file):
         else:
             raise FileNotFoundError('%s input file not found'% file)
         return indata
+
+    @property
+    def savefile(self):
+        simfolder = self.get_folder_name()
+        return os.path.join(simfolder, 'job.npy')
+
     
 def read_inputfile(file):
     """
@@ -793,7 +832,8 @@ def read_inputfile(file):
             #create calculations
             for combo in combos:
                 calc = Calculation.generate(indata)
-                calc.add_from_dict(ci, keys=["mode", "pair_style", "pair_coeff", "pair_style_options", "npt", "repeat", "n_equilibration_steps",
+                calc.add_from_dict(ci, keys=["mode", "pair_style", "pair_coeff", "pair_style_options", "npt", 
+                                "repeat", "n_equilibration_steps",
                                 "n_switching_steps", "n_print_steps", "n_iterations", "potential_file", "spring_constants",
                                 "melting_cycle", "equilibration_control", "folder_prefix", "temperature_high"])
                 calc.lattice = combo[0]["lattice"]
@@ -803,3 +843,13 @@ def read_inputfile(file):
                 calc.temperature = combo[2]
                 calculations.append(calc)
     return calculations
+
+
+def save_job(job):
+    filename = os.path.join(job.simfolder, 'job.npy')
+    np.save(filename, job)
+
+def load_job(filename):
+    job = np.load(filename, allow_pickle=True).flatten()[0]
+    return job
+
diff --git a/calphy/kernel.py b/calphy/kernel.py
index 5197ba1..dac1083 100644
--- a/calphy/kernel.py
+++ b/calphy/kernel.py
@@ -55,29 +55,84 @@ def run_jobs(inputfile):
     calculations = read_inputfile(inputfile)
     print("Total number of %d calculations found" % len(calculations))
 
-    for count, calc in enumerate(calculations):
-
-        identistring = calc.create_identifier()
-        scriptpath = os.path.join(os.getcwd(), ".".join([identistring, "sub"]))
-        errfile = os.path.join(os.getcwd(), ".".join([identistring, "err"]))
-
-        #the below part assigns the schedulers
-        #now we have to write the submission scripts for the job
-        #parse Queue and import module
-        if calc.queue.scheduler == "local":
-            scheduler = pq.Local(calc.queue.__dict__, cores=calc.queue.cores)
-        elif calc.queue.scheduler == "slurm":
-            scheduler = pq.SLURM(calc.queue.__dict__, cores=calc.queue.cores)
-        elif calc.queue.scheduler == "sge":
-            scheduler = pq.SGE(calc.queue.__dict__, cores=calc.queue.cores)
-        else:
-            raise ValueError("Unknown scheduler")
-
-        #for lattice just provide the number of position
-        scheduler.maincommand = "calphy_kernel -i %s -k %d"%(inputfile, 
-            count)
-        scheduler.write_script(scriptpath)
-        _ = scheduler.submit()
+    if calculations[0].script_mode:
+        for count, calc in enumerate(calculations):
+            identistring = calc.create_identifier()
+            scriptpath = os.path.join(os.getcwd(), ".".join([identistring, "sub"]))
+            errfile = os.path.join(os.getcwd(), ".".join([identistring, "err"]))
+            if calc.queue.scheduler == "local":
+                scheduler = pq.Local(calc.queue.__dict__, cores=calc.queue.cores)
+            elif calc.queue.scheduler == "slurm":
+                scheduler = pq.SLURM(calc.queue.__dict__, cores=calc.queue.cores)
+            elif calc.queue.scheduler == "sge":
+                scheduler = pq.SGE(calc.queue.__dict__, cores=calc.queue.cores)
+            else:
+                raise ValueError("Unknown scheduler")
+
+            #add run commands
+            command = f'calphy_run_averaging -i {inputfile} -k {count}'
+            scheduler.queueoptions['commands'].append(command)
+
+            command = f'cp input.yaml {os.path.join(os.getcwd(), identistring)}'
+            scheduler.queueoptions['commands'].append(command)
+
+            command = f'cd {os.path.join(os.getcwd(), identistring)}'
+            scheduler.queueoptions['commands'].append(command)
+            if calc.queue.cores > 1:
+                #here turn on mpi
+                command = f'{calc.mpi_executable} -np {calc.queue.cores} {calc.lammps_executable} -in averaging.lmp'
+            else:
+                command = f'{calc.lammps_executable} -in averaging.lmp'
+            scheduler.queueoptions['commands'].append(command)
+            scheduler.queueoptions['commands'].append('cd ..')
+
+            command = f'calphy_process_averaging -i {inputfile} -k {count}'
+            scheduler.queueoptions['commands'].append(command)
+
+            command = f'calphy_run_integration -i {inputfile} -k {count}'
+            scheduler.queueoptions['commands'].append(command)
+
+            command = f'cd {os.path.join(os.getcwd(), identistring)}'
+            scheduler.queueoptions['commands'].append(command)
+            if calc.queue.cores > 1:
+                #here turn on mpi
+                command = f'{calc.mpi_executable} -np {calc.queue.cores} {calc.lammps_executable} -in integration.lmp'
+            else:
+                command = f'{calc.lammps_executable} -in integration.lmp'
+            scheduler.queueoptions['commands'].append(command)
+            scheduler.queueoptions['commands'].append('cd ..')
+
+            command = f'calphy_process_integration -i {inputfile} -k {count}'
+            scheduler.maincommand = command
+            scheduler.write_script(scriptpath)
+            #_ = scheduler.submit()
+
+
+
+    else:
+        for count, calc in enumerate(calculations):
+
+            identistring = calc.create_identifier()
+            scriptpath = os.path.join(os.getcwd(), ".".join([identistring, "sub"]))
+            errfile = os.path.join(os.getcwd(), ".".join([identistring, "err"]))
+
+            #the below part assigns the schedulers
+            #now we have to write the submission scripts for the job
+            #parse Queue and import module
+            if calc.queue.scheduler == "local":
+                scheduler = pq.Local(calc.queue.__dict__, cores=calc.queue.cores)
+            elif calc.queue.scheduler == "slurm":
+                scheduler = pq.SLURM(calc.queue.__dict__, cores=calc.queue.cores)
+            elif calc.queue.scheduler == "sge":
+                scheduler = pq.SGE(calc.queue.__dict__, cores=calc.queue.cores)
+            else:
+                raise ValueError("Unknown scheduler")
+
+            #for lattice just provide the number of position
+            scheduler.maincommand = "calphy_kernel -i %s -k %d"%(inputfile, 
+                count)
+            scheduler.write_script(scriptpath)
+            _ = scheduler.submit()
 
 
 
diff --git a/calphy/lattice.py b/calphy/lattice.py
index 02285ec..2a4f7d7 100644
--- a/calphy/lattice.py
+++ b/calphy/lattice.py
@@ -26,6 +26,7 @@
 from pylammpsmpi import LammpsLibrary
 import numpy as np
 import pyscal.core as pc
+from ase.io import read
 
 """
 Conversion factors for creating initial lattices
@@ -98,30 +99,35 @@ def check_dump_file(infile):
     except:
         return None
 
-def check_data_file(infile):
+def check_data_file(infile, script_mode=False):
     try:
-        lmp = LammpsLibrary(cores=1, 
-            working_directory=os.getcwd())
-        lmp.units("metal")
-        lmp.boundary("p p p")
-        lmp.atom_style("atomic")
-        lmp.timestep(0.001)            
-        lmp.read_data(infile)
-        natoms = lmp.natoms
-        #now we convert to a dump file and read the concentration
-        trajfile = ".".join([infile, "dump"])
-        lmp.command("mass * 1.0")
-        lmp.dump("2 all custom", 1, trajfile,"id type x y z")
-        lmp.run(0)
-        lmp.undump(2)
+        if not script_mode:
+            lmp = LammpsLibrary(cores=1, 
+                working_directory=os.getcwd())
+            lmp.units("metal")
+            lmp.boundary("p p p")
+            lmp.atom_style("atomic")
+            lmp.timestep(0.001)            
+            lmp.read_data(infile)
+            natoms = lmp.natoms
+            #now we convert to a dump file and read the concentration
+            trajfile = ".".join([infile, "dump"])
+            lmp.command("mass * 1.0")
+            lmp.dump("2 all custom", 1, trajfile,"id type x y z")
+            lmp.run(0)
+            lmp.undump(2)
+            lmp.close()
+            format = 'lammps-dump'
+        else:
+            trajfile = read(infile, format='lammps-data', style='atomic')
+            format = 'ase'
         #now use pyscal to read it in,
         sys = pc.System()
-        sys.read_inputfile(trajfile)
+        sys.read_inputfile(trajfile, format=format)
         atoms = sys.atoms
         types = [atom.type for atom in atoms]
         xx, xxcounts = np.unique(types, return_counts=True)
         conc = xxcounts/np.sum(xxcounts)
-        lmp.close()
         return (natoms, conc)
     except:
         return None
diff --git a/calphy/liquid.py b/calphy/liquid.py
index 77f02dd..db9c358 100644
--- a/calphy/liquid.py
+++ b/calphy/liquid.py
@@ -145,11 +145,9 @@ def run_averaging(self):
         self.dump_current_snapshot(lmp, "traj.equilibration_stage1.dat")
         self.check_if_solidfied(lmp, "traj.equilibration_stage1.dat")
         self.dump_current_snapshot(lmp, "traj.equilibration_stage2.dat")
-        lmp = ph.write_data(lmp, "current.data")
+        lmp = ph.write_data(lmp, "conf.equilibration.data")
         lmp.close()
 
-        #process the trajectory
-        self.process_traj("traj.equilibration_stage2.dat", "conf.equilibration.data")
 
 
 
diff --git a/calphy/phase.py b/calphy/phase.py
index 7b5f7ed..6523cb0 100644
--- a/calphy/phase.py
+++ b/calphy/phase.py
@@ -50,7 +50,7 @@ class Phase:
     """
     def __init__(self, calculation=None, simfolder=None):
 
-        self.calc = calculation
+        self.calc = copy.deepcopy(calculation)
         self.simfolder = simfolder
         
         logfile = os.path.join(self.simfolder, "calphy.log")
@@ -164,6 +164,7 @@ def __init__(self, calculation=None, simfolder=None):
         self.ly = None
         self.lz = None
 
+
         #now manually tune pair styles
         if self.calc.pair_style is not None:
             self.logger.info("pair_style: %s"%self.calc.pair_style_with_options[0])
@@ -447,6 +448,28 @@ def run_pressure_convergence(self, lmp):
         -------
         None
 
+        Notes
+        -----
+        Take the equilibrated structure and rigorously check for pressure convergence.
+        The cycle is stopped when the average pressure is within the given cutoff of the target pressure.
+        """
+        if self.calc.script_mode:
+            self.run_minimal_pressure_convergence(lmp)
+        else:
+            self.run_iterative_pressure_convergence(lmp)
+
+    def run_iterative_pressure_convergence(self, lmp):
+        """
+        Run a pressure convergence routine
+
+        Parameters
+        ----------
+        lmp: LAMMPS object
+
+        Returns
+        -------
+        None
+
         Notes
         -----
         Take the equilibrated structure and rigorously check for pressure convergence.
@@ -507,8 +530,52 @@ def run_pressure_convergence(self, lmp):
 
         lmp.command("unfix            2")
 
+    def run_minimal_pressure_convergence(self, lmp):
+        """
+        Run a pressure convergence routine
+
+        Parameters
+        ----------
+        lmp: LAMMPS object
+
+        Returns
+        -------
+        None
+
+        Notes
+        -----
+        Take the equilibrated structure and rigorously check for pressure convergence.
+        The cycle is stopped when the average pressure is within the given cutoff of the target pressure.
+        """
+
+        #apply fixes
+        if self.calc.equilibration_control == "nose-hoover":
+            self.fix_nose_hoover(lmp)
+        else:
+            self.fix_berendsen(lmp)
+
+
+        lmp.command("fix              2 all ave/time %d %d %d v_mlx v_mly v_mlz v_mpress file avg.dat"%(int(self.calc.md.n_every_steps),
+            int(self.calc.md.n_repeat_steps), int(self.calc.md.n_every_steps*self.calc.md.n_repeat_steps)))
+        
+        lmp.command("run              %d"%int(self.calc.md.n_small_steps))
+        lmp.command("run              %d"%int(self.calc.n_equilibration_steps))
+
+        #unfix thermostat and barostat
+        if self.calc.equilibration_control == "nose-hoover":
+            self.unfix_nose_hoover(lmp)
+        else:
+            self.unfix_berendsen(lmp)
+
+        lmp.command("unfix            2")
 
     def run_constrained_pressure_convergence(self, lmp):
+        if self.calc.script_mode:
+            self.run_minimal_constrained_pressure_convergence(lmp)
+        else:
+            self.run_iterative_constrained_pressure_convergence(lmp)
+
+    def run_iterative_constrained_pressure_convergence(self, lmp):
         """
         """
         lmp.command("velocity         all create %f %d"%(self.calc._temperature, np.random.randint(0, 10000)))
@@ -524,34 +591,19 @@ def run_constrained_pressure_convergence(self, lmp):
         converged = False
         for i in range(int(self.calc.md.n_cycles)):
             lmp.command("run              %d"%int(self.calc.md.n_small_steps))
-            ncount = int(self.calc.md.n_small_steps)//int(self.calc.md.n_every_steps*self.calc.md.n_repeat_steps)
-            #now we can check if it converted
-            file = os.path.join(self.simfolder, "avg.dat")
-            lx, ly, lz, ipress = np.loadtxt(file, usecols=(1, 2, 3, 4), unpack=True)
             
-            lxpc = ipress
-            mean = np.mean(lxpc)
-            std = np.std(lxpc)
-            volatom = np.mean((lx*ly*lz)/self.natoms)
+            #now we can check if it converted
+            mean, std,volatom = self.process_pressure()
             self.logger.info("At count %d mean pressure is %f with %f vol/atom"%(i+1, mean, volatom))
 
             if (np.abs(mean - lastmean)) < 50*self.calc.tolerance.pressure:
                 #here we actually have to set the pressure
-                self.calc._pressure = mean
-                std = np.std(lxpc)
-                volatom = np.mean((lx*ly*lz)/self.natoms)
-                self.logger.info("At count %d mean pressure is %f with %f vol/atom"%(i+1, mean, volatom))
-                self.lx = np.round(np.mean(lx[-ncount+1:]), decimals=3)
-                self.ly = np.round(np.mean(ly[-ncount+1:]), decimals=3)
-                self.lz = np.round(np.mean(lz[-ncount+1:]), decimals=3)
-                self.volatom = volatom
-                self.vol = self.lx*self.ly*self.lz
-                self.logger.info("finalized vol/atom %f at pressure %f"%(self.volatom, mean))
-                self.logger.info("Avg box dimensions x: %f, y: %f, z:%f"%(self.lx, self.ly, self.lz))
-                #now run for msd
+                self.finalise_pressure()
                 converged = True
                 break
+
             lastmean = mean
+        
         lmp.command("unfix            1")
         lmp.command("unfix            2")
 
@@ -559,11 +611,69 @@ def run_constrained_pressure_convergence(self, lmp):
             lmp.close()
             raise ValueError("pressure did not converge")
 
+    def process_pressure(self,):
+        if self.calc.script_mode:
+            ncount = int(self.calc.n_equilibration_steps)//int(self.calc.md.n_every_steps*self.calc.md.n_repeat_steps)
+        else: 
+            ncount = int(self.calc.md.n_small_steps)//int(self.calc.md.n_every_steps*self.calc.md.n_repeat_steps)
+        
+        #now we can check if it converted
+        file = os.path.join(self.simfolder, "avg.dat")
+        lx, ly, lz, ipress = np.loadtxt(file, usecols=(1, 2, 3, 4), unpack=True)
+        lxpc = ipress
+        mean = np.mean(lxpc)
+        std = np.std(lxpc)
+        volatom = np.mean((lx*ly*lz)/self.natoms)
+        return mean, std, volatom
+
+    def finalise_pressure(self,):
+        if self.calc.script_mode:
+            ncount = int(self.calc.n_equilibration_steps)//int(self.calc.md.n_every_steps*self.calc.md.n_repeat_steps)
+        else: 
+            ncount = int(self.calc.md.n_small_steps)//int(self.calc.md.n_every_steps*self.calc.md.n_repeat_steps)
+
+        file = os.path.join(self.simfolder, "avg.dat")
+        lx, ly, lz, ipress = np.loadtxt(file, usecols=(1, 2, 3, 4), unpack=True)        
+        lxpc = ipress
+        mean = np.mean(lxpc)
+        std = np.std(lxpc)
+        volatom = np.mean((lx*ly*lz)/self.natoms)
+
+        self.calc._pressure = mean
+        self.lx = np.round(np.mean(lx[-ncount+1:]), decimals=3)
+        self.ly = np.round(np.mean(ly[-ncount+1:]), decimals=3)
+        self.lz = np.round(np.mean(lz[-ncount+1:]), decimals=3)
+        self.volatom = volatom
+        self.vol = self.lx*self.ly*self.lz
+        self.logger.info("finalized vol/atom %f at pressure %f"%(self.volatom, mean))
+        self.logger.info("Avg box dimensions x: %f, y: %f, z:%f"%(self.lx, self.ly, self.lz))
+
+
+    def run_minimal_constrained_pressure_convergence(self, lmp):
+        """
+        """
+        lmp.command("velocity         all create %f %d"%(self.calc._temperature, np.random.randint(0, 10000)))
+        lmp.command("fix              1 all nvt temp %f %f %f"%(self.calc._temperature, self.calc._temperature, self.calc.md.thermostat_damping[1]))
+        lmp.command("thermo_style     custom step pe press vol etotal temp lx ly lz")
+        lmp.command("thermo           10")
+        
+        #this is when the averaging routine starts
+        lmp.command("fix              2 all ave/time %d %d %d v_mlx v_mly v_mlz v_mpress file avg.dat"%(int(self.calc.md.n_every_steps),
+            int(self.calc.md.n_repeat_steps), int(self.calc.md.n_every_steps*self.calc.md.n_repeat_steps)))
+
+        lmp.command("run              %d"%int(self.calc.md.n_small_steps))
+        lmp.command("run              %d"%int(self.calc.n_equilibration_steps))
+
+        lmp.command("unfix            1")
+        lmp.command("unfix            2")
+
     def process_traj(self, filename, outfilename):
         """
         Process the out trajectory after averaging cycle and 
         extract a configuration to run integration
 
+        Will be phased out...
+
         Parameters
         ----------
         None
@@ -577,7 +687,9 @@ def process_traj(self, filename, outfilename):
         files = ptp.split_trajectory(trajfile)
         conf = os.path.join(self.simfolder, outfilename)
 
-        ph.reset_timestep(conf, os.path.join(self.simfolder, "current.data"), init_commands=self.calc.md.init_commands)
+        ph.reset_timestep(conf, os.path.join(self.simfolder, "current.data"), 
+            init_commands=self.calc.md.init_commands,
+            script_mode=self.calc.script_mode)
 
         os.remove(trajfile)
         for file in files:
@@ -758,11 +870,12 @@ def reversible_scaling(self, iteration=1):
         lmp.command("run               %d"%self.calc.n_equilibration_steps)
 
         #check melting or freezing
-        self.dump_current_snapshot(lmp, "traj.temp.dat")
-        if solid:
-            self.check_if_melted(lmp, "traj.temp.dat")
-        else:
-            self.check_if_solidfied(lmp, "traj.temp.dat")
+        if not self.calc.script_mode:
+            self.dump_current_snapshot(lmp, "traj.temp.dat")
+            if solid:
+                self.check_if_melted(lmp, "traj.temp.dat")
+            else:
+                self.check_if_solidfied(lmp, "traj.temp.dat")
 
         lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
 
@@ -897,11 +1010,12 @@ def temperature_scaling(self, iteration=1):
         lmp.command("run               0")
         lmp.command("undump            2")
         
-        self.dump_current_snapshot(lmp, "traj.temp.dat")
-        if solid:
-            self.check_if_melted(lmp, "traj.temp.dat")
-        else:
-            self.check_if_solidfied(lmp, "traj.temp.dat")
+        if not self.calc.script_mode:
+            self.dump_current_snapshot(lmp, "traj.temp.dat")
+            if solid:
+                self.check_if_melted(lmp, "traj.temp.dat")
+            else:
+                self.check_if_solidfied(lmp, "traj.temp.dat")
 
         #start reverse loop
         lmp.command("variable          lambda equal ramp(${lf},${li})")
@@ -1017,4 +1131,4 @@ def integrate_pressure_scaling(self, return_values=False):
             nsims=self.calc.n_iterations, return_values=return_values)
 
         if return_values:
-            return res
\ No newline at end of file
+            return res    
\ No newline at end of file
diff --git a/calphy/queuekernel.py b/calphy/queuekernel.py
index 4b887eb..f1e19c1 100644
--- a/calphy/queuekernel.py
+++ b/calphy/queuekernel.py
@@ -122,13 +122,7 @@ def main():
     calc = calculations[kernel]
     
     #format and parse the arguments
-    identistring = calc.create_identifier()
-    simfolder = os.path.join(os.getcwd(), identistring)
-
-    #if folder exists, delete it -> then create
-    if os.path.exists(simfolder):
-        shutil.rmtree(simfolder)
-    os.mkdir(simfolder)
+    simfolder = calc.create_folders()
 
     if calc.mode == "melting_temperature":
         os.rmdir(simfolder)
diff --git a/calphy/solid.py b/calphy/solid.py
index f402302..c2636bf 100644
--- a/calphy/solid.py
+++ b/calphy/solid.py
@@ -57,6 +57,12 @@ def __init__(self, calculation=None, simfolder=None):
         simfolder=simfolder)
 
     def run_spring_constant_convergence(self, lmp):
+        if self.calc.script_mode:
+            self.run_minimal_spring_constant_convergence(lmp)
+        else:
+            self.run_iterative_spring_constant_convergence(lmp)
+
+    def run_iterative_spring_constant_convergence(self, lmp):
         """
         """
         lmp.command("fix              3 all nvt temp %f %f %f"%(self.calc._temperature, self.calc._temperature, self.calc.md.thermostat_damping[1]))
@@ -69,38 +75,13 @@ def run_spring_constant_convergence(self, lmp):
             laststd = 0.00
             for i in range(self.calc.md.n_cycles):
                 lmp.command("run              %d"%int(self.calc.md.n_small_steps))
-                ncount = int(self.calc.md.n_small_steps)//int(self.calc.md.n_every_steps*self.calc.md.n_repeat_steps)
-                #now we can check if it converted
-                file = os.path.join(self.simfolder, "msd.dat")
-                quant = np.loadtxt(file, usecols=(1,), unpack=True)[-ncount+1:]
-                quant = 3*kb*self.calc._temperature/quant
-                #self.logger.info(quant)
-                mean = np.mean(quant)
-                std = np.std(quant)
-                self.logger.info("At count %d mean k is %f std is %f"%(i+1, mean, std))
-                if (np.abs(laststd - std) < self.calc.tolerance.spring_constant):
+                k_mean, k_std = self.analyse_spring_constants()
+                self.logger.info("At count %d mean k is %f std is %f"%(i+1, k_mean[0], k_std[0]))
+                if (np.abs(laststd - k_std[0]) < self.calc.tolerance.spring_constant):
                     #now reevaluate spring constants
-                    k = []
-                    for i in range(self.calc.n_elements):
-                        quant = np.loadtxt(file, usecols=(i+1, ), unpack=True)[-ncount+1:]
-                        quant = 3*kb*self.calc._temperature/quant
-                        k.append(np.round(np.mean(quant), decimals=2))
-
-                    #first replace any provided values with user values
-                    if ph.check_if_any_is_not_none(self.calc.spring_constants):
-                        spring_constants = copy.copy(self.calc.spring_constants)
-                        k = ph.replace_nones(spring_constants, k, logger=self.logger)
-                    
-                    #add sanity checks
-                    k = ph.validate_spring_constants(k, logger=self.logger)
-                    
-                    #now save
-                    self.k = k
-
-                    self.logger.info("finalized sprint constants")
-                    self.logger.info(self.k)
+                    self.assign_spring_constants(k_mean)                    
                     break
-                laststd = std
+                laststd = k_std
 
         else:
             if not (len(self.calc.spring_constants) == self.calc.n_elements):
@@ -114,6 +95,68 @@ def run_spring_constant_convergence(self, lmp):
 
         lmp.command("unfix         3")
 
+    def analyse_spring_constants(self):
+        """
+        Analyse spring constant routine
+        """
+        if self.calc.script_mode:
+            ncount = int(self.calc.md.n_small_steps)//int(self.calc.md.n_every_steps*self.calc.md.n_repeat_steps)
+        else:
+            ncount = int(self.calc.n_equilibration_steps)//int(self.calc.md.n_every_steps*self.calc.md.n_repeat_steps)
+        
+        #now we can check if it converted
+        file = os.path.join(self.simfolder, "msd.dat")
+        k_mean = []
+        k_std = []
+        for i in range(self.calc.n_elements):
+            quant = np.loadtxt(file, usecols=(i+1, ), unpack=True)[-ncount+1:]
+            quant = 3*kb*self.calc._temperature/quant
+            k_mean.append(np.round(np.mean(quant), decimals=2))
+            k_std.append(np.round(np.std(quant), decimals=2))
+        return k_mean, k_std
+        
+
+    def assign_spring_constants(self, k):
+        """
+        Here the spring constants are finalised, add added to the class
+        """
+        #first replace any provided values with user values
+        if ph.check_if_any_is_not_none(self.calc.spring_constants):
+            spring_constants = copy.copy(self.calc.spring_constants)
+            k = ph.replace_nones(spring_constants, k, logger=self.logger)
+        
+        #add sanity checks
+        k = ph.validate_spring_constants(k, logger=self.logger)
+        
+        #now save
+        self.k = k
+
+        self.logger.info("finalized sprint constants")
+        self.logger.info(self.k)
+
+
+    def run_minimal_spring_constant_convergence(self, lmp):
+        """
+        """
+        lmp.command("fix              3 all nvt temp %f %f %f"%(self.calc._temperature, self.calc._temperature, self.calc.md.thermostat_damping[1]))
+        
+        #apply fix
+        lmp = ph.compute_msd(lmp, self.calc)
+        
+        if ph.check_if_any_is_none(self.calc.spring_constants):
+            #similar averaging routine
+                lmp.command("run              %d"%int(self.calc.md.n_small_steps))
+                lmp.command("run              %d"%int(self.calc.n_equilibration_steps))
+
+        else:
+            if not (len(self.calc.spring_constants) == self.calc.n_elements):
+                raise ValueError("Spring constant input length should be same as number of elements, spring constant length %d, # elements %d"%(len(self.calc.spring_constants), self.calc.n_elements))
+
+            #still run a small NVT cycle
+            lmp.command("run              %d"%int(self.calc.md.n_small_steps))
+
+        lmp.command("unfix         3")
+
     def run_averaging(self):
         """
         Run averaging routine
@@ -136,8 +179,38 @@ def run_averaging(self):
         is calculated.
         At the end of the run, the averaged box dimensions are calculated. 
         """
+        if self.calc.script_mode:
+            self.run_minimal_averaging()
+        else:
+            self.run_interactive_averaging()
+
+    def run_interactive_averaging(self):
+        """
+        Run averaging routine
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+
+        Notes
+        -----
+        Run averaging routine using LAMMPS. Starting from the initial lattice two different routines can
+        be followed:
+        If pressure is specified, MD simulations are run until the pressure converges within the given
+        threshold value.
+        If `fix_lattice` option is True, then the input structure is used as it is and the corresponding pressure
+        is calculated.
+        At the end of the run, the averaged box dimensions are calculated. 
+        """
+
         lmp = ph.create_object(self.cores, self.simfolder, self.calc.md.timestep, 
-            self.calc.md.cmdargs, self.calc.md.init_commands)
+            self.calc.md.cmdargs, 
+            init_commands=self.calc.md.init_commands,
+            script_mode=self.calc.script_mode)
 
         #set up structure
         lmp = ph.create_structure(lmp, self.calc, species=self.calc.n_elements+self.calc._ghost_element_count)
@@ -182,12 +255,90 @@ def run_averaging(self):
         #check for melting
         self.dump_current_snapshot(lmp, "traj.equilibration_stage2.dat")
         self.check_if_melted(lmp, "traj.equilibration_stage2.dat")
-        lmp = ph.write_data(lmp, "current.data")
+        lmp = ph.write_data(lmp, "conf.equilibration.data")
         #close object and process traj
         lmp.close()
-        self.process_traj("traj.equilibration_stage2.dat", "conf.equilibration.data")
 
 
+    def run_minimal_averaging(self):
+        """
+        Run averaging routine
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+
+        Notes
+        -----
+        Run averaging routine using LAMMPS. Starting from the initial lattice two different routines can
+        be followed:
+        If pressure is specified, MD simulations are run until the pressure converges within the given
+        threshold value.
+        If `fix_lattice` option is True, then the input structure is used as it is and the corresponding pressure
+        is calculated.
+        At the end of the run, the averaged box dimensions are calculated. 
+        """
+        lmp = ph.create_object(self.cores, self.simfolder, self.calc.md.timestep, 
+            self.calc.md.cmdargs, 
+            init_commands=self.calc.md.init_commands,
+            script_mode=self.calc.script_mode)
+
+        #set up structure
+        lmp = ph.create_structure(lmp, self.calc, species=self.calc.n_elements+self.calc._ghost_element_count)
+
+        #set up potential
+        if self.calc.potential_file is None:
+            lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
+        else:
+            lmp.command("include %s"%self.calc.potential_file)
+
+        #add some computes
+        lmp.command("variable         mvol equal vol")
+        lmp.command("variable         mlx equal lx")
+        lmp.command("variable         mly equal ly")
+        lmp.command("variable         mlz equal lz")
+        lmp.command("variable         mpress equal press")
+
+        #Run if a constrained lattice is not needed
+        if not self.calc._fix_lattice:
+            if self.calc._pressure == 0:
+                self.run_zero_pressure_equilibration(lmp)
+            else:
+                self.run_finite_pressure_equilibration(lmp)
+
+
+            #this is when the averaging routine starts
+            self.run_pressure_convergence(lmp)
+
+            #dump snapshot and check if melted
+            self.dump_current_snapshot(lmp, "traj.equilibration_stage1.dat")
+        
+        #run if a constrained lattice is used
+        else:
+            #routine in which lattice constant will not varied, but is set to a given fixed value
+            self.run_constrained_pressure_convergence(lmp)
+
+        #start MSD calculation routine
+        #there two possibilities here - if spring constants are provided, use it. If not, calculate it
+        self.run_spring_constant_convergence(lmp)
+
+        #check for melting
+        self.dump_current_snapshot(lmp, "traj.equilibration_stage2.dat")
+        lmp = ph.write_data(lmp, "conf.equilibration.data")
+        #close object and process traj
+
+        #now serialise script
+        file = os.path.join(self.simfolder, 'averaging.lmp')
+        lmp.write(file)
+
+    def process_averaging_results(self):
+        k_m, k_s = self.analyse_spring_constants()
+        self.assign_spring_constants(k_m)
+        self.finalise_pressure()
 
     def run_integration(self, iteration=1):
         """
@@ -208,7 +359,9 @@ def run_integration(self, iteration=1):
         the lambda parameter. See algorithm 4 in publication.
         """
         lmp = ph.create_object(self.cores, self.simfolder, self.calc.md.timestep, 
-            self.calc.md.cmdargs, self.calc.md.init_commands)
+            self.calc.md.cmdargs, 
+            init_commands=self.calc.md.init_commands,
+            script_mode=self.calc.script_mode)
 
         #read in the conf file
         #conf = os.path.join(self.simfolder, "conf.equilibration.dump")
@@ -307,7 +460,11 @@ def run_integration(self, iteration=1):
         lmp.command("unfix             f4")
 
         #close object
-        lmp.close()
+        if not self.calc.script_mode:
+            lmp.close()
+        else:
+            file = os.path.join(self.simfolder, 'integration.lmp')
+            lmp.write(file)
 
 
     def thermodynamic_integration(self):
diff --git a/examples/example_01/input.yaml b/examples/example_01/input.yaml
index 3121260..aaf873c 100644
--- a/examples/example_01/input.yaml
+++ b/examples/example_01/input.yaml
@@ -18,10 +18,10 @@ md:
   timestep: 0.001
   thermostat_damping: 0.1
   barostat_damping: 0.1
-  init_commands:
-     - timestep 0.002
-     - neighbor 0.6 bin
-  
+  equilibration_control: nose_hoover
+
+tolerance:
+  pressure: 1  
 queue:
   scheduler: local
   cores: 4
diff --git a/notebooks/running_from_notebook.ipynb b/notebooks/running_from_notebook.ipynb
index 561110a..dc2b186 100644
--- a/notebooks/running_from_notebook.ipynb
+++ b/notebooks/running_from_notebook.ipynb
@@ -306,7 +306,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.2"
+   "version": "3.10.9"
   }
  },
  "nbformat": 4,
diff --git a/setup.py b/setup.py
index ebf6cba..201a90a 100644
--- a/setup.py
+++ b/setup.py
@@ -63,6 +63,10 @@
         'console_scripts': [
             'calphy = calphy.kernel:main',
             'calphy_kernel = calphy.queuekernel:main',
+            'calphy_run_averaging = calphy.clitools:run_averaging',
+            'calphy_process_averaging = calphy.clitools:process_averaging',
+            'calphy_run_integration = calphy.clitools:run_integration',
+            'calphy_process_integration = calphy.clitools:process_integration',
         ],
     }
 )