diff --git a/sbmy_control/parameters_card.py b/sbmy_control/parameters_card.py
index bbdc1a6..5384adb 100644
--- a/sbmy_control/parameters_card.py
+++ b/sbmy_control/parameters_card.py
@@ -47,6 +47,15 @@ def register_arguments_card(parser:ArgumentParser):
     parser.add_argument("--OutputFCsDensity", type=str, default=None, help="Output FCs density file.")
     parser.add_argument("--OutputFCsSnapshot", type=str, default=None, help="Output FCs snapshot file.")
     parser.add_argument("--OutputRngStateLPT", type=str, default=None, help="Output RNG state file.")
+    ## Tests with phiBCs and density
+    parser.add_argument("--WriteGravPot", type=bool, default=False, help="Write gravitational potential.")
+    parser.add_argument("--OutputGravitationalPotentialBase", type=str, default=None, help="Output gravitational potential base.")
+    parser.add_argument("--MeshGravPot", type=int, default=None, help="Mesh for gravitational potential.")
+    parser.add_argument("--WriteDensity", type=bool, default=False, help="Write density.")
+    parser.add_argument("--OutputDensityBase", type=str, default=None, help="Output density base.")
+    parser.add_argument("--MeshDensity", type=int, default=None, help="Mesh for density.")
+    parser.add_argument("--LoadPhiBCs", type=bool, default=False, help="Load phiBCs.")
+    parser.add_argument("--InputPhiBCsBase", type=str, default=None, help="Input phiBCs file base.")
     parser.add_argument("--WriteReferenceFrame", type=bool, default=False, help="Write reference frame (COCA).")
     parser.add_argument("--OutputMomentaBase", type=str, default=None, help="Output momenta base (COCA).")
     parser.add_argument("--ReadReferenceFrame", type=bool, default=False, help="Read reference frame (COCA).")
@@ -126,6 +135,16 @@ def parse_arguments_card(parsed_args):
             OutputFCsDensity=parsed_args.OutputFCsDensity,
             OutputFCsSnapshot=parsed_args.OutputFCsSnapshot,
             OutputRngStateLPT=parsed_args.OutputRngStateLPT,
+            ## Tests with phiBCs and density
+            WriteGravPot=parsed_args.WriteGravPot,
+            OutputGravitationalPotentialBase=parsed_args.OutputGravitationalPotentialBase,
+            MeshGravPot=parsed_args.MeshGravPot,
+            WriteDensity=parsed_args.WriteDensity,
+            OutputDensityBase=parsed_args.OutputDensityBase,
+            MeshDensity=parsed_args.MeshDensity,
+            LoadPhiBCs=parsed_args.LoadPhiBCs,
+            InputPhiBCsBase=parsed_args.InputPhiBCsBase,
+            ## Cosmological parameters
             h=cosmo_dict["h"],
             Omega_m=cosmo_dict["Omega_m"],
             Omega_b=cosmo_dict["Omega_b"],
@@ -211,6 +230,17 @@ def parse_arguments_card(parsed_args):
         card_dict["OutputFCsSnapshot"] = main_dict["resultdir"]+"final_particles_"+main_dict["simname"]+".gadget3"
     if card_dict["OutputRngStateLPT"] is None:
         card_dict["OutputRngStateLPT"] = main_dict["workdir"]+"rng_state.h5"
+    ## Tests with phiBCs and density
+    if card_dict["OutputGravitationalPotentialBase"] is None:
+        card_dict["OutputGravitationalPotentialBase"] = main_dict["workdir"]+"gravpot_"+main_dict["simname"]
+    if card_dict["MeshGravPot"] is None:
+        card_dict["MeshGravPot"] = card_dict["N_PM_mesh"]
+    if card_dict["OutputDensityBase"] is None:
+        card_dict["OutputDensityBase"] = main_dict["workdir"]+"density_"+main_dict["simname"]
+    if card_dict["MeshDensity"] is None:
+        card_dict["MeshDensity"] = card_dict["N_PM_mesh"]
+    if card_dict["InputPhiBCsBase"] is None:
+        card_dict["InputPhiBCsBase"] = main_dict["workdir"]+"gravpot_tCOLA"
     if card_dict["OutputMomentaBase"] is None:
         card_dict["OutputMomentaBase"] = main_dict["workdir"]+"momenta_"+main_dict["simname"]+"_"
     if card_dict["InputMomentaBase"] is None:
@@ -333,6 +363,15 @@ def create_parameter_card_dict(
                     OutputFCsDensity:str = 'fcs_density.h5',
                     OutputFCsSnapshot:str = 'fcs_particles.gadget3',
 
+                    ## Tests with phiBCs and density
+                    WriteGravPot:bool = True,
+                    OutputGravitationalPotentialBase:str = 'gravitational_potential.h5',
+                    MeshGravPot:int = 128,
+                    WriteDensity:bool = False,
+                    OutputDensityBase:str = 'density.h5',
+                    MeshDensity:int = 128,
+                    LoadPhiBCs:bool = False,
+                    InputPhiBCsBase:str = 'gravitational_potential.h5',
                     ## COCA parameters
                     WriteReferenceFrame:bool = False,
                     OutputMomentaBase:str = 'momenta_',
@@ -403,6 +442,17 @@ def create_parameter_card_dict(
             OutputLPTPotential1=OutputLPTPotential1,
             OutputLPTPotential2=OutputLPTPotential2,
             OutputTilesBase=OutputTilesBase,
+
+            # Tests with phiBCs and density
+            WriteGravPot=int(WriteGravPot),
+            OutputGravitationalPotentialBase=OutputGravitationalPotentialBase,
+            MeshGravPot=MeshGravPot,
+            WriteDensity=int(WriteDensity),
+            OutputDensityBase=OutputDensityBase,
+            MeshDensity=MeshDensity,
+            LoadPhiBCs=int(LoadPhiBCs),
+            InputPhiBCsBase=InputPhiBCsBase,
+
             h=h,
             Omega_m=Omega_m,
             Omega_b=Omega_b,