diff --git a/src/simulation/Simulation.cpp b/src/simulation/Simulation.cpp
index 82fce6b7e..be970ee2b 100644
--- a/src/simulation/Simulation.cpp
+++ b/src/simulation/Simulation.cpp
@@ -31,6 +31,52 @@ void Simulation::Load(const GameSave *save, bool includePressure, Vec2<int> bloc
 
 	RecalcFreeParticles(false);
 
+	struct ExistingParticle
+	{
+		int id;
+		Vec2<int> pos;
+	};
+	std::vector<ExistingParticle> existingParticles;
+	auto pasteArea = RES.OriginRect() & RectSized(partP, save->blockSize * CELL);
+	for (int i = 0; i <= parts_lastActiveIndex; i++)
+	{
+		if (parts[i].type)
+		{
+			auto p = Vec2<int>{ int(parts[i].x + 0.5f), int(parts[i].y + 0.5f) };
+			if (pasteArea.Contains(p))
+			{
+				existingParticles.push_back({ i, p });
+			}
+		}
+	}
+	std::sort(existingParticles.begin(), existingParticles.end(), [](const auto &lhs, const auto &rhs) {
+		return std::tie(lhs.pos.Y, lhs.pos.X) < std::tie(rhs.pos.Y, rhs.pos.X);
+	});
+	PlaneAdapter<std::vector<size_t>> existingParticleIndices(pasteArea.Size(), existingParticles.size());
+	{
+		auto lastPos = Vec2<int>{ -1, -1 }; // not a valid pos in existingParticles
+		for (auto it = existingParticles.begin(); it != existingParticles.end(); ++it)
+		{
+			if (lastPos != it->pos)
+			{
+				existingParticleIndices[it->pos - pasteArea.TopLeft] = it - existingParticles.begin();
+				lastPos = it->pos;
+			}
+		}
+	}
+	auto removeExistingParticles = [this, pasteArea, &existingParticles, &existingParticleIndices](Vec2<int> p) {
+		auto rp = p - pasteArea.TopLeft;
+		if (existingParticleIndices.Size().OriginRect().Contains(rp))
+		{
+			auto index = existingParticleIndices[rp];
+			for (auto it = existingParticles.begin() + index; it != existingParticles.end() && it->pos == p; ++it)
+			{
+				kill_part(it->id);
+			}
+			existingParticleIndices[rp] = existingParticles.size();
+		}
+	};
+
 	std::map<unsigned int, unsigned int> soapList;
 	for (int n = 0; n < NPART && n < save->particlesCount; n++)
 	{
@@ -69,9 +115,6 @@ void Simulation::Load(const GameSave *save, bool includePressure, Vec2<int> bloc
 			continue;
 		}
 
-		// Mark location to be cleaned of existing particles.
-		pmap[y][x] = -1;
-
 		if (elements[tempPart.type].CreateAllowed)
 		{
 			if (!(*(elements[tempPart.type].CreateAllowed))(this, -3, int(tempPart.x + 0.5f), int(tempPart.y + 0.5f), tempPart.type))
@@ -80,6 +123,8 @@ void Simulation::Load(const GameSave *save, bool includePressure, Vec2<int> bloc
 			}
 		}
 
+		removeExistingParticles({ x, y });
+
 		// Allocate particle (this location is guaranteed to be empty due to "full scan" logic above)
 		if (pfree == -1)
 			break;
@@ -171,35 +216,12 @@ void Simulation::Load(const GameSave *save, bool includePressure, Vec2<int> bloc
 		{
 			parts[i].tmp3 = 0;
 		}
-
-		if (!parts[i].type)
-		{
-			continue;
-		}
-
-		// Mark to be preserved in the loop below.
-		parts[i].type |= 1 << PMAPBITS;
 	}
 	parts_lastActiveIndex = NPART-1;
 	force_stacking_check = true;
 	Element_PPIP_ppip_changed = 1;
 
-	// Loop through particles to find particles in need of being killed
-	for (int i = 0; i <= parts_lastActiveIndex; i++)
-	{
-		if (parts[i].type)
-		{
-			int x = int(parts[i].x + 0.5f);
-			int y = int(parts[i].y + 0.5f);
-			bool preserve = parts[i].type & (1 << PMAPBITS);
-			parts[i].type &= ~(1 << PMAPBITS);
-			if (pmap[y][x] == -1 && !preserve)
-			{
-				kill_part(i);
-			}
-		}
-	}
-
+	// Sort out pmap, just to be on the safe side.
 	RecalcFreeParticles(false);
 
 	// fix SOAP links using soapList, a map of old particle ID -> new particle ID