In recent years, conventional ultrasound (US) imaging devices have been adapted with the photoacoustic (PA) imaging capabilities to simultaneously provide both anatomical and molecular optical contrasts of soft biological tissues. To help optimize the design parameters of such dual modality imaging devices, we present a numerical simulation approach for Bmode beamformed US and multispectral PA imaging using a linear ultrasound transducer array surrounded by a light source. We combined the finite element based simulation platforms for ultrasound and light propagation, K-wave and NIRFast respectively, to model the ultrasound and photoacoustic effects in deep tissue, and created an effective hybrid platform for simulating US and multispectral PA imaging of different configurations. We also developed and applied a spectral unmixing algorithm on multispectral photoacoustic images, obtained from multiple optical wavelengths, to map different molecules (e.g., Indocyanogreen (ICG), Deoxyhemoglobin (Hb), and Oxyhemoglobin (HbO2)) present inside the tissue background. The multi-spectral plots and unmixed spectral images clearly delineated the molecular contrast arising from different regions inside the tissue. The presented simulation platform allows for optimization of key design parameters of both US and PA imaging devices, such as the size of ultrasonic transducer array, and size and the distribution of light sources. Our results demonstrate that the ability to mimic the imaging performance of such dual modality deep tissue-imaging device will help to achieve high molecular sensitivity for the targeted clinical application, thus functioning as a powerful tool for medical device design.