To write a 4-spin term, you won’t be able to use the square lattice helper function (the one that gives you b.s1 and b.s2), since that function gives you an array only of nearest-neighbor pairs of sites, not 4-site groups.
Therefore, you will need to write your own logic that loops over the bonds and 4-site groups (I assume plaquettes?) of the lattice you want to do, mapping the 2D coordinates to the 1D site indices of the MPO and MPS that is going into your DMRG calculation. I would suggest printing out the b.s1 and b.s2 values for the nearest-neighbor case to study what that does, and draw figures showing the 2D to 1D mapping, and looking at the other posts on this forum asking about similar questions to study how this is done. For example this discussion.