The accelerated expansion of the Universe can be explained by either the cosmological constant Λ, a repulsive dark energy, or modifications to General Relativity (GR) on cosmological scales. The key to distinguishing these possibilities is to study the growth of large-scale structures over cosmic time, measured primarily by three methods: cosmic shear from weak lensing (WL), redshift-distortion (RSD) of galaxy clustering, and the traditional cluster cosmology using the abundance of galaxy clusters as function of mass. Although cosmic shear and RSD are commonly studied using galaxies, the two effects are much more enhanced around galaxy clusters. In this dissertation, I develop novel yet powerful techniques to combine cluster WL, galaxy infall kinematics, and cluster abundances to derive constraints on cosmological models and test of modified gravity (MG) theories.One of the main systematical uncertainties of cluster cosmology is in cluster mass calibration. Common mass observables, including galaxy richness, X-ray luminosity, and Sunyaev-Zel'dovich (SZ) decrements, correlate with cluster true mass with some scatter. Using MaxBCG cluster catalog derived from the Sloan Digital Sky Survey (SDSS) and assuming Λ-Cold Dark Matter (ΛCDM) model, I developed an alternative method for constraining cosmological parameters, i.e., the matter density Ωm and the clustering amplitude ς8, while simultaneously calibrating the scatter in the richness-mass relation, by combining the abundance of clusters as function of richness and the large-scale cluster WL profiles. We find ς8(Ωm/0.325)0.501=0.828±0.049, consistent with constraints from other MaxBCG studies that use WL measurements on small scales. The (Ωm,ς8) constraint is consistent with and orthogonal to the one inferred from WMAP CMB data, reflecting agreement with the structure growth predicted by General Relativity for a ΛCDM cosmological model. A joint constraint assuming ΛCDM yields Ωm=0.298-0.020+0.019 and ς8=0.831-0.020+0.020. For these parameters and our best-fit scatter we obtain a tightly constrained mean richness-mass relation of MaxBCG clusters, N200=25.4(M/3.61 ×1014 h-1Msun)0.74, with a normalization uncertainty of 1.5%.Cluster WL is currently the only viable method to measure cluster mass profiles beyond the virial radius. I developed an novel method to infer the dynamical cluster mass profiles by extracting the characteristic infall velocities as function of radius, from measurements of the redshift-space cluster-galaxy cross-correlation function, ξcgs. Using the Millennium simulation, I calibrate the galaxy kinematic profiles around clusters as a two-component mixture model comprised of a virialized component with an isotropic velocity distribution and a infall component described by a skewed t-distribution with separate radial and tangential dispersions. By convolving with the real-space cluster-galaxy cross-correlation function ξcgs, I show that the galaxy infall kinematics (GIK) model provides accurate predictions for ξcgs, and I show that the measurement of ξcgs can be inverted to recover the GIK profiles for any given cluster sample. In particular, the inferred characteristic infall velocity profile is a promising tool for estimating the dynamical mass of clusters. As a first-cut experiment, I apply the GIK-modelling to ξcgs measured for rich galaxy groups derived from the SDSS main galaxy sample, and the inferred masses agree well with those estimated from virial scaling and weak lensing.To explain cosmic acceleration, most MG theories rely on an extra scalar field on top of the metric tensor field, mediating a "fifth" force on intermediate scales. However, photons do not respond to the scalar field, so there will be a mismatch between lensing and dynamical mass estimates if gravity is modified on relevant scales. Since GR has to be recovered in high density regions to evade stringent tests in the solar system and binary pulsars, MG theories generally invoke some "screening" mechanisms to either make the effective mass of the scalar field environment-dependent (Chameleon mechanism) or suppress the scalar field close to sources (Vainshtein mechanism). I investigate the impact of MG on the GIK profiles around massive clusters using two suites of MG simulations, namely, the f(R) and Galileon models, which employ the Chameleon and Vainshtein screening mechanisms, respectively. For clusters of M~1014h-1Msun at z~0.25, I find ~100-200 km/s enhancement in the characteristic infall velocity and ~50-100 km/s broadening in the radial and tangential velocity dispersions at r~5 Mpc compared to GR predictions. These deviations from GR are detectable via the GIK-modeling of ξcgs measured from existing and future galaxy redshift surveys.For future work, I will address the major source of systematic uncertainties in GIK calibration, the imperfect understanding of the impact of galaxy formation physics on GIK, by comparing GIK profiles measured from mock galaxy catalogs constructed using different Halo Occupation Distribution (HOD) prescriptions and semi-analytic models (SAMs) within large-volume, high-resolution cosmological simulations. By calibrating the GIK as function of halo mass in simulations of different cosmological parameters, I will be able to combine WL, GIK, and abundances of galaxy clusters as one of the most powerful diagnostics of theories of cosmic acceleration and modified gravity.