While I am aware that VASP and QE are widely used for solid-state materials and PySCF/OpenMolcas are primarily designed for molecular systems.
I want to know why should I use PySCF/OpenMolcas for molecular calculations over VASP/QE and whether those packages could replace VASP/QE for periodic systems?
The main difference is that PySCF, OpenMolcas, and other popular software such as Gaussian, ORCA, Molpro, Psi4, MRCC, CFOUR, Dalton, GAMESS, COLUMBUS, ACES, etc., are very well-developed for using Gaussian basis sets. In a typical molecule, there will be nuclei at certain irregular places, and the electronic wavefunction will look very different near those nuclei compared to away from those nuclei, and Gaussian basis sets that are centered at each nucleus, offer an excellent way to help model the shape of the wavefunction when we have such irregularities at very specific places.
In PySCF and OpenMOLCAS, we also have access to well-developed and extremely high-accuracy (compared to VASP and QuantumESPRESSO) wavefunction-based methods for solving the Schroedinger equation: full-configuration-interaction (FCI), coupled-cluster, casscf and in OpenMolcas we even have methods built on top of CASSCF such as CASPT2, RASPT2, and GASPT2. FCI offers an exact solution to the Schroedinger equation for a given basis set, so not only is the Gaussian basis set capturing the fine details of the shape of the equation more accurately than a plane-wave basis set typically would, the actual solution to the Schroedinger equation based on that basis set is also capable to be more accurate in the Gaussian-based software.
VASP and Quantum ESPRESSO use plane-wave basis sets which don't capture the local irregularities of the wavefunction near the nuclei as well as Gaussian basis sets do. Furthermore, the capabilities to do calculations such as FCI, coupled cluster, etc. are not nearly as well-developed in VASP and Quantum ESPRESSO. VASP and Quantum ESPRESSO are designed to model comparatively huge materials with far more nuclei and electrons, so the methods available are mainly low-cost methods such as DFT with plane-wave basis sets: these will not be able to compete in accuracy with the above-mentioned wavefunction-based methods with Gaussian basis sets, and when modeling smaller molecules, the standards of accuracy tend to be a lot higher (it's not okay to just get a band gap accurate to within ±0.1 eV of the exact solution to the Schroedinger equation, the goal will be much closer to ±1 cm-1 accuracy or at least ±1 kcal/mol accuracy). A plane-wave basis set can certainly compete with a Guassian basis sets for accuracy if the former contains far more functions than the latter, but you would need orders of magnitude more functions, so the vast majority of the people that model molecules do not use plane-wave basis sets.
PySCF has some functionality to do calculations with periodic boundary conditions, but it is not nearly as developed as the functionality for that in VASP and QuantumESPRESSO. Contrarily, I do not yet recommend using VASP or Quantum ESPRESSO for calculations with Gaussian basis sets.
