An empirical method has been designed to account for the van der Waals interactions in practical molecular calculations with density functional theory. For each atom pair separated at a distance R, the method adds to the density functional electronic structure calculations an additional attraction energy E-vdW=-f(d)(R)C6R-6, where f(d)(R) is the damping function which equals to one at large value of R and zero at small value of R. The coefficients C-6 for pair interactions between hydrogen, carbon, nitrogen, and oxygen atoms have been developed in this work by a least-square fitting to the molecular C-6 coefficients obtained from the dipole oscillator strength distribution method by Meath and co-workers. Two forms of the damping functions have been studied, with one dropping to zero at short distances much faster than the other. Four density functionals have been examined: Becke's three parameter hybrid functional with the Lee-Yang-Parr correlation functional, Becke's 1988 exchange functional with the LYP correlation functional, Becke's 1988 exchange functional with Perdew and Wang's 1991 (PW91) correlation functional, and PW91 exchange and correlation functional. The method has been applied to three systems where the van der Waals attractions are known to be important: rare-gas diatomic molecules, stacking of base pairs and polyalanines' conformation stabilities. The results show that this empirical method, with the damping function dropping to zero smoothly, provides a significant correction to both of the Becke's hybrid functional and the PW91 exchange and correlation functional. Results are comparable to the corresponding second-order Moller-Plesset calculations in many cases.