diff --git a/Docs/source/dataanalysis/inputs.2d b/Docs/source/dataanalysis/inputs.2d
index 389766f13e6..4c5fe7efb58 100644
--- a/Docs/source/dataanalysis/inputs.2d
+++ b/Docs/source/dataanalysis/inputs.2d
@@ -13,7 +13,8 @@ insitu.int = 2
 insitu.start = 0
 insitu.config = ez2d.xml
 geometry.coord_sys   = 0                  # 0: Cartesian
-geometry.is_periodic = 0 0      # Is periodic?
+boundary.field_lo = pml pml
+boundary.field_hi = pml pml
 geometry.prob_lo     = -125.e-6 -149.e-6 
 geometry.prob_hi     =  125.e-6    1.e-6 
 warpx.fine_tag_lo    =  -12.e-6 -110.e-6
@@ -29,7 +30,6 @@ algo.particle_pusher = 0
 algo.particle_shape = 3
 algo.maxwell_solver = ckc
 warpx.use_filter = 1
-warpx.do_pml = 1
 warpx.pml_ncell = 10
 warpx.cfl = .99
 warpx.do_moving_window = 1
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index f1bf915ba1c..6f13a972825 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -211,7 +211,7 @@ Domain Boundary Conditions
 --------------------------
 
 * ``boundary.field_lo`` and ``boundary.field_hi`` (`2 strings` for 2D, `3 strings` for 3D, `pml` by default)
-    Boundary conditions applied to field at the lower and upper domain boundaries. Depending on the type of boundary condition, the value for ``geometry.is_periodic`` will be set, overriding the user-input for the input parameter, ``geometry.is_periodic``. If not set, the default value for the fields at the domain boundary will be set to pml.
+    Boundary conditions applied to fields at the lower and upper domain boundaries.
     Options are:
 
     * ``Periodic``: This option can be used to set periodic domain boundaries. Note that if the fields for lo in a certain dimension are set to periodic, then the corresponding upper boundary must also be set to periodic. If particle boundaries are not specified in the input file, then particles boundaries by default will be set to periodic. If particles boundaries are specified, then they must be set to periodic corresponding to the periodic field boundaries.
diff --git a/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_mr.py b/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_mr.py
index 8ea29032da9..4dc144b3eb6 100644
--- a/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_mr.py
+++ b/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_mr.py
@@ -35,7 +35,6 @@
                         hi = [25e-6, 25e-6, 200.e-6])
 
 solver = picmi.ElectromagneticSolver(grid=grid, cfl=1,
-                                     warpx_do_pml = True,
                                      warpx_pml_ncell = 10)
 
 beam_distribution = picmi.UniformDistribution(density = 1.e23,
diff --git a/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rz_multimode_analyze.py b/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rz_multimode_analyze.py
index 16a3f05752a..22ab1679b24 100644
--- a/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rz_multimode_analyze.py
+++ b/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rz_multimode_analyze.py
@@ -91,7 +91,7 @@
                              moving_window_zvelocity = 0.,
                              warpx_max_grid_size=64)
 
-solver = picmi.ElectromagneticSolver(grid=grid, cfl=1., warpx_do_pml=False)
+solver = picmi.ElectromagneticSolver(grid=grid, cfl=1.)
 
 ##########################
 # diagnostics
diff --git a/Examples/Tests/SilverMueller/inputs_2d_x b/Examples/Tests/SilverMueller/inputs_2d_x
index bbfa1aad16e..12251d037c8 100644
--- a/Examples/Tests/SilverMueller/inputs_2d_x
+++ b/Examples/Tests/SilverMueller/inputs_2d_x
@@ -23,7 +23,6 @@ warpx.verbose = 1
 
 # Algorithms
 warpx.cfl = 1.0
-warpx.do_pml = 0
 warpx.use_filter = 0
 
 warpx.do_moving_window = 0
diff --git a/Examples/Tests/SilverMueller/inputs_2d_z b/Examples/Tests/SilverMueller/inputs_2d_z
index cf3ff31471b..ef6b5d01f48 100644
--- a/Examples/Tests/SilverMueller/inputs_2d_z
+++ b/Examples/Tests/SilverMueller/inputs_2d_z
@@ -24,7 +24,6 @@ warpx.verbose = 1
 
 # Algorithms
 warpx.cfl = 1.0
-warpx.do_pml = 0
 warpx.use_filter = 0
 
 warpx.do_moving_window = 0
diff --git a/Examples/Tests/SilverMueller/inputs_rz_z b/Examples/Tests/SilverMueller/inputs_rz_z
index d3308dac74e..228fb8a1180 100644
--- a/Examples/Tests/SilverMueller/inputs_rz_z
+++ b/Examples/Tests/SilverMueller/inputs_rz_z
@@ -27,7 +27,6 @@ warpx.n_rz_azimuthal_modes = 2
 
 # Algorithms
 warpx.cfl = 1.0
-warpx.do_pml = 0
 warpx.use_filter = 0
 
 warpx.do_moving_window = 0
diff --git a/Examples/Tests/boundaries/inputs_3d b/Examples/Tests/boundaries/inputs_3d
index 749ce51a567..2fba4a6d016 100644
--- a/Examples/Tests/boundaries/inputs_3d
+++ b/Examples/Tests/boundaries/inputs_3d
@@ -22,8 +22,6 @@ boundary.field_hi = pec pec periodic
 boundary.particle_lo = reflecting absorbing periodic
 boundary.particle_hi = reflecting absorbing periodic
 
-warpx.do_pml = 0
-
 # Algorithms
 algo.particle_shape = 1
 
diff --git a/Examples/Tests/galilean/inputs_rz b/Examples/Tests/galilean/inputs_rz
index 2704e00fe86..92f84a3ab25 100644
--- a/Examples/Tests/galilean/inputs_rz
+++ b/Examples/Tests/galilean/inputs_rz
@@ -36,7 +36,6 @@ particles.species_names = electrons ions
 
 warpx.do_nodal = 1
 warpx.use_filter = 1
-warpx.do_pml = 0
 
 psatd.nox = 16
 psatd.noy = 16
diff --git a/Examples/Tests/multi_J/inputs_rz b/Examples/Tests/multi_J/inputs_rz
index 1248a174a6f..588d907fbdf 100644
--- a/Examples/Tests/multi_J/inputs_rz
+++ b/Examples/Tests/multi_J/inputs_rz
@@ -23,7 +23,6 @@ algo.particle_shape = 3
 
 # Numerics
 warpx.do_moving_window = 1
-warpx.do_pml = 0
 warpx.moving_window_dir = z
 warpx.moving_window_v = 1.
 warpx.n_rz_azimuthal_modes = 1
diff --git a/Examples/Tests/photon_pusher/analysis_photon_pusher.py b/Examples/Tests/photon_pusher/analysis_photon_pusher.py
index 2347d0baca6..83dfb0c6a8a 100755
--- a/Examples/Tests/photon_pusher/analysis_photon_pusher.py
+++ b/Examples/Tests/photon_pusher/analysis_photon_pusher.py
@@ -120,10 +120,10 @@ def generate():
         f.write("amr.max_grid_size = 8\n")
         f.write("amr.plot_int = 1\n")
         f.write("geometry.coord_sys   = 0\n")
-        f.write("geometry.is_periodic = 1 1 1\n")
+        f.write("boundary.field_lo = periodic periodic periodic\n")
+        f.write("boundary.field_hi = periodic periodic periodic\n")
         f.write("geometry.prob_lo = -0.5e-6 -0.5e-6 -0.5e-6\n")
         f.write("geometry.prob_hi = 0.5e-6 0.5e-6 0.5e-6\n")
-        f.write("warpx.do_pml = 0\n")
         f.write("algo.charge_deposition = standard\n")
         f.write("algo.field_gathering = energy-conserving\n")
         f.write("warpx.cfl = 1.0\n")
diff --git a/Examples/Tests/plasma_lens/inputs_3d b/Examples/Tests/plasma_lens/inputs_3d
index c2a96a7a1b7..67b549eb29d 100644
--- a/Examples/Tests/plasma_lens/inputs_3d
+++ b/Examples/Tests/plasma_lens/inputs_3d
@@ -16,7 +16,6 @@ boundary.field_hi = pec pec pec
 boundary.particle_lo = absorbing absorbing absorbing
 boundary.particle_hi = absorbing absorbing absorbing
 
-warpx.do_pml = 0
 warpx.const_dt = 1.e-7
 warpx.do_electrostatic = labframe
 
diff --git a/Examples/Tests/radiation_reaction/test_const_B_analytical/analysis_classicalRR.py b/Examples/Tests/radiation_reaction/test_const_B_analytical/analysis_classicalRR.py
index c008527a9ed..a09654b5372 100755
--- a/Examples/Tests/radiation_reaction/test_const_B_analytical/analysis_classicalRR.py
+++ b/Examples/Tests/radiation_reaction/test_const_B_analytical/analysis_classicalRR.py
@@ -163,10 +163,10 @@ def generate():
         f.write("amr.blocking_factor = 32\n")
         f.write("amr.max_grid_size = 64\n")
         f.write("geometry.coord_sys   = 0\n")
-        f.write("geometry.is_periodic = 1 1 1\n")
+        f.write("boundary.field_lo = periodic periodic periodic\n")
+        f.write("boundary.field_hi = periodic periodic periodic\n")
         f.write("geometry.prob_lo = {} {} {}\n".format(-sim_size, -sim_size, -sim_size))
         f.write("geometry.prob_hi = {} {} {}\n".format(sim_size, sim_size, sim_size))
-        f.write("warpx.do_pml = 0\n")
         f.write("algo.charge_deposition = standard\n")
         f.write("algo.field_gathering = energy-conserving\n")
         f.write("warpx.cfl = 1.0\n")
diff --git a/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs_3d b/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs_3d
index 65a2eb4b8a2..5b78921f1bf 100644
--- a/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs_3d
+++ b/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs_3d
@@ -10,7 +10,6 @@ amr.max_grid_size = 64
 geometry.coord_sys   = 0
 geometry.prob_lo = -8e-07 -8e-07 -8e-07
 geometry.prob_hi = 8e-07 8e-07 8e-07
-warpx.do_pml = 0
 algo.charge_deposition = standard
 algo.field_gathering = energy-conserving
 warpx.cfl = 1.0
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
index f9495d2b551..dff586261ae 100644
--- a/Python/pywarpx/picmi.py
+++ b/Python/pywarpx/picmi.py
@@ -525,7 +525,6 @@ class ElectromagneticSolver(picmistandard.PICMI_ElectromagneticSolver):
     def init(self, kw):
         assert self.method is None or self.method in ['Yee', 'CKC', 'PSATD'], Exception("Only 'Yee', 'CKC', and 'PSATD' are supported")
 
-        self.do_pml = kw.pop('warpx_do_pml', None)
         self.pml_ncell = kw.pop('warpx_pml_ncell', None)
 
         if self.method == 'PSATD':
@@ -539,7 +538,6 @@ def initialize_inputs(self):
 
         self.grid.initialize_inputs()
 
-        pywarpx.warpx.do_pml = self.do_pml
         pywarpx.warpx.pml_ncell = self.pml_ncell
         pywarpx.warpx.do_nodal = self.l_nodal
 
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 1880bcc311c..311116edc55 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -433,6 +433,7 @@ public:
     void DampJPML (int lev, PatchType patch_type);
 
     void CopyJPML ();
+    bool isAnyBoundaryPML();
 
     /**
      * \brief Synchronize the nodal points of the PML MultiFabs
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index ab66cdcc187..d31382c2022 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -666,14 +666,13 @@ WarpX::ReadParameters ()
             }
         }
 
-        pp_warpx.query("do_pml", do_pml);
         pp_warpx.query("pml_ncell", pml_ncell);
         pp_warpx.query("pml_delta", pml_delta);
         pp_warpx.query("pml_has_particles", pml_has_particles);
         pp_warpx.query("do_pml_j_damping", do_pml_j_damping);
         pp_warpx.query("do_pml_in_domain", do_pml_in_domain);
 
-        if (do_multi_J && do_pml)
+        if (do_multi_J && isAnyBoundaryPML())
         {
             amrex::Abort("Multi-J algorithm not implemented with PMLs");
         }
@@ -737,38 +736,8 @@ WarpX::ReadParameters ()
         }
 
 #ifdef WARPX_DIM_RZ
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( do_pml==0,
-            "PML are not implemented in RZ geometry ; please set `warpx.do_pml=0`");
-#endif
-        // setting default to 0
-        Vector<int> parse_do_pml_Lo(AMREX_SPACEDIM,0);
-        // Switching pml lo to 1 when do_pml = 1 and if domain is non-periodic
-        // Note to remove this code when new BC API is fully functional
-        if (do_pml == 1) {
-            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-                if ( Geom(0).isPeriodic(idim) == 0) parse_do_pml_Lo[idim] = 1;
-            }
-        }
-        pp_warpx.queryarr("do_pml_Lo", parse_do_pml_Lo);
-        do_pml_Lo[0] = parse_do_pml_Lo[0];
-        do_pml_Lo[1] = parse_do_pml_Lo[1];
-#if (AMREX_SPACEDIM == 3)
-        do_pml_Lo[2] = parse_do_pml_Lo[2];
-#endif
-        // setting default to 0
-        Vector<int> parse_do_pml_Hi(AMREX_SPACEDIM,0);
-        // Switching pml hi to 1 when do_pml = 1 and if domain is non-periodic
-        // Note to remove this code when new BC API is fully functional
-        if (do_pml == 1) {
-            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-                if ( Geom(0).isPeriodic(idim) == 0) parse_do_pml_Hi[idim] = 1;
-            }
-        }
-        pp_warpx.queryarr("do_pml_Hi", parse_do_pml_Hi);
-        do_pml_Hi[0] = parse_do_pml_Hi[0];
-        do_pml_Hi[1] = parse_do_pml_Hi[1];
-#if (AMREX_SPACEDIM == 3)
-        do_pml_Hi[2] = parse_do_pml_Hi[2];
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( isAnyBoundaryPML() == false,
+            "PML are not implemented in RZ geometry; please set a different boundary condition using boundary.field_lo and boundary.field_hi.");
 #endif
 
         if ( (do_pml_j_damping==1)&&(do_pml_in_domain==0) ){
@@ -1234,7 +1203,9 @@ WarpX::BackwardCompatibility ()
         amrex::Abort("warpx.use_kspace_filter is not supported anymore. "
                      "Please use the flag use_filter, see documentation.");
     }
-
+    if ( pp_warpx.query("do_pml", backward_int) ) {
+        amrex::Abort( "do_pml is not supported anymore. Please use boundary.field_lo and boundary.field_hi to set the boundary conditions.");
+    }
     ParmParse pp_interpolation("interpolation");
     if (pp_interpolation.query("nox", backward_int) ||
         pp_interpolation.query("noy", backward_int) ||
@@ -2305,3 +2276,13 @@ WarpX::PicsarVersion ()
     return std::string("Unknown");
 #endif
 }
+
+bool
+WarpX::isAnyBoundaryPML()
+{
+    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+        if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::PML) return true;
+        if ( WarpX::field_boundary_hi[idim] == FieldBoundaryType::PML) return true;
+    }
+    return false;
+}
diff --git a/Tools/PerformanceTests/automated_test_1_uniform_rest_32ppc b/Tools/PerformanceTests/automated_test_1_uniform_rest_32ppc
index ba38b055f23..1a31a9aaefe 100644
--- a/Tools/PerformanceTests/automated_test_1_uniform_rest_32ppc
+++ b/Tools/PerformanceTests/automated_test_1_uniform_rest_32ppc
@@ -23,7 +23,6 @@ boundary.particle_hi = absorbing absorbing periodic
 warpx.verbose = 1
 
 algo.particle_shape = 3
-warpx.do_pml = 0
 
 # CFL
 warpx.cfl = 1.0
diff --git a/Tools/PerformanceTests/automated_test_2_uniform_rest_1ppc b/Tools/PerformanceTests/automated_test_2_uniform_rest_1ppc
index 7ff6eb6ddce..17b7a02b43f 100644
--- a/Tools/PerformanceTests/automated_test_2_uniform_rest_1ppc
+++ b/Tools/PerformanceTests/automated_test_2_uniform_rest_1ppc
@@ -23,7 +23,6 @@ boundary.particle_hi = absorbing absorbing absorbing
 warpx.verbose = 1
 
 algo.particle_shape = 3
-warpx.do_pml = 1
 
 # CFL
 warpx.cfl = 1.0
diff --git a/Tools/PerformanceTests/automated_test_3_uniform_drift_4ppc b/Tools/PerformanceTests/automated_test_3_uniform_drift_4ppc
index 107cac09bdf..7a671167f6c 100644
--- a/Tools/PerformanceTests/automated_test_3_uniform_drift_4ppc
+++ b/Tools/PerformanceTests/automated_test_3_uniform_drift_4ppc
@@ -24,7 +24,6 @@ warpx.verbose = 1
 
 # Algorithms
 algo.particle_shape = 3
-warpx.do_pml = 0
 
 # CFL
 warpx.cfl = 1.0
diff --git a/Tools/PerformanceTests/automated_test_4_labdiags_2ppc b/Tools/PerformanceTests/automated_test_4_labdiags_2ppc
index 705880a039a..f8337c8bae7 100644
--- a/Tools/PerformanceTests/automated_test_4_labdiags_2ppc
+++ b/Tools/PerformanceTests/automated_test_4_labdiags_2ppc
@@ -24,7 +24,6 @@ warpx.verbose = 1
 algo.particle_shape = 3
 warpx.use_filter = 1
 warpx.cfl = 1.0
-warpx.do_pml = 0
 
 # Moving window
 warpx.do_moving_window = 1
diff --git a/Tools/PerformanceTests/automated_test_5_loadimbalance b/Tools/PerformanceTests/automated_test_5_loadimbalance
index f004588033d..bdfd874ef69 100644
--- a/Tools/PerformanceTests/automated_test_5_loadimbalance
+++ b/Tools/PerformanceTests/automated_test_5_loadimbalance
@@ -22,7 +22,6 @@ boundary.particle_hi = absorbing absorbing periodic
 warpx.verbose = 1
 algo.load_balance_intervals = -5
 algo.particle_shape = 3
-warpx.do_pml = 0
 
 # CFL
 warpx.cfl = 1.0
diff --git a/Tools/PerformanceTests/automated_test_6_output_2ppc b/Tools/PerformanceTests/automated_test_6_output_2ppc
index 36c2cb01c91..89d69c476da 100644
--- a/Tools/PerformanceTests/automated_test_6_output_2ppc
+++ b/Tools/PerformanceTests/automated_test_6_output_2ppc
@@ -23,7 +23,6 @@ boundary.particle_hi = absorbing absorbing absorbing
 warpx.verbose = 1
 
 algo.particle_shape = 3
-warpx.do_pml = 0
 
 # CFL
 warpx.cfl = 1.0