code: def solvate_bilayer(self):
          
          solvent_number_density = self.solvent.n_particles / np.prod(self.solvent.periodicity)
  
          lengths = self.lipid_box.lengths
          water_box_z = self.solvent_per_layer / (lengths[0] * lengths[1] * solvent_number_density)
  
          mins = self.lipid_box.mins
          maxs = self.lipid_box.maxs
          bilayer_solvent_box = mb.Box(mins=[mins[0], mins[1], maxs[2]],
                                       maxs=[maxs[0], maxs[1], maxs[2] + 2 * water_box_z])
  
          self.solvent_components.add(mb.fill_box(self.solvent, bilayer_solvent_box))
docs: Solvate the constructed bilayer.