diff --git a/calphy/input.py b/calphy/input.py
index d053518..d9bb622 100644
--- a/calphy/input.py
+++ b/calphy/input.py
@@ -61,6 +61,7 @@ def __init__(self):
         self._fix_lattice = False
         self._pair_style = None
         self._pair_coeff = None
+        self._potential_file = None
         self._reference_phase = None
         self._lattice_constant = 0
         self._repeat = [1, 1, 1]
@@ -316,6 +317,17 @@ def pair_coeff(self, val):
         val = self.fix_paths(val)
         self._pair_coeff = val
 
+    @property
+    def potential_file(self):
+        return self._potential_file
+
+    @potential_file.setter
+    def potential_file(self, val):
+        if os.path.exists(val):
+            self._potential_file = val
+        else:
+            raise FileNotFoundError("File %s not found"%val)
+
     @property
     def reference_phase(self):
         return self._reference_phase
@@ -577,7 +589,7 @@ def read_inputfile(file):
             for combo in combos:
                 calc = Calculation.generate(indata)
                 calc.from_dict(ci, keys=["mode", "pair_style", "pair_coeff", "repeat", "n_equilibration_steps",
-                                "n_switching_steps", "n_print_steps", "n_iterations"])
+                                "n_switching_steps", "n_print_steps", "n_iterations", "potential_file"])
                 calc.lattice = combo[0]["lattice"]
                 calc.lattice_constant = combo[0]["lattice_constant"]
                 calc.reference_phase = combo[0]["reference_phase"]
diff --git a/calphy/solid.py b/calphy/solid.py
index ead2cc5..3d191a6 100644
--- a/calphy/solid.py
+++ b/calphy/solid.py
@@ -82,7 +82,10 @@ def run_averaging(self):
         lmp = ph.create_structure(lmp, self.calc)
 
         #set up potential
-        lmp = ph.set_potential(lmp, self.calc)
+        if self.calc.potential_file is None:
+            lmp = ph.set_potential(lmp, self.calc)
+        else:
+            lmp.command("include %s"%self.calc.potential_file)
 
         #add some computes
         lmp.command("variable         mvol equal vol")
@@ -310,7 +313,10 @@ def run_integration(self, iteration=1):
         lmp = ph.read_dump(lmp, conf, species=self.calc.n_elements)
 
         #set up potential
-        lmp = ph.set_potential(lmp, self.calc)
+        if self.calc.potential_file is None:
+            lmp = ph.set_potential(lmp, self.calc)
+        else:
+            lmp.command("include %s"%self.calc.potential_file)
 
         #remap the box to get the correct pressure
         lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)
diff --git a/docs/source/documentation/inputfile.md b/docs/source/documentation/inputfile.md
index a0245f8..00db3e8 100644
--- a/docs/source/documentation/inputfile.md
+++ b/docs/source/documentation/inputfile.md
@@ -11,7 +11,7 @@ The inputfile is `yaml` formatted. In this section the possible keys in the inpu
 | :-: | :-: | :-: | :-: | :-: |
 | [mode](#mode) | [lattice](#lattice) | [reference_phase](#reference_phase) | [temperature](#temperature) | [pressure](#pressure) |
 | [temperature_high](#temperature_high) | [lattice_constant](#lattice_constant) | [repeat](#repeat) | [n_iterations](#n_iterations) | [n_switching_steps](#n_switching_steps) | 
-| [n_equilibration_steps](#n_equilibration_steps) | [pair_style](#pair_style) | [pair_coeff](#pair_coeff) | [n_print_steps](#n_print_steps) | 
+| [n_equilibration_steps](#n_equilibration_steps) | [pair_style](#pair_style) | [pair_coeff](#pair_coeff) | [n_print_steps](#n_print_steps) | [potential_file](#potential_file) | 
 
 | `md` block  | | | | |
 | :-: | :-: | :-: | :-: | :-: |
@@ -279,6 +279,19 @@ The [pair coeff](https://lammps.sandia.gov/doc/pair_coeff.html) command for LAMM
 
 ---
 
+#### <a name="potential_file"></a>`potential_file`
+
+_type_: string      
+_default_: None  
+_example_:
+```
+pair_coeff: "/home/calc/potential.inp"
+```
+
+If specified, the `pair_style` and `pair_coeff` commands are not used, but rather the potential is read in from the provided input file using `include` command in LAMMPS. This allows the use of more complex or multiple potential files. Due to the `hybrid/scaled` styles employed in calphy, **this option only works with mode `fe` and `reference_phase` solid.**
+
+---
+
 #### <a name="n_equilibration_steps"></a>`n_equilibration_steps`
 
 _type_: int