The Lattice Boltzmann (LB) method has evolved from a theoretical novelty to a ubiquitous, versatile and powerful computational methodology for both fundamental research and engineering applications. It is a kinetic-based mesoscopic approach that bridges the microscales and macroscales, which offers distinctive advantages in simulation fidelity and computational efficiency. Since the emergence of the LB method, its application in multiphase flows has always been a very important theme of the method. With the development in the past two decades, many multiphase LB models have been proposed. After a review of potential approaches, the pseudopotential method introduced by Shan and Chen is implemented in an in-house fully-differentiable PyTorch-based Lattice Boltzmann simulation framework. To overcome the shortcomings of the original forcing scheme, a robust scheme based on Guo's approach is implemented, resulting in a significant improvement in the thermodynamic consistency of the method. Moreover, the pseudopotential method is extended to incorporate realistic equations of state such as the Carnahan-Starling EOS. Finally, potential performance bottlenecks are identified by profiling the simulation framework.