We perform molecular dynamics simulations of screening by bound target electrons in low energy nuclear reactions. Quantum effects corresponding to the Pauli and Heisenberg principle are enforced by constraints. We show that the enhancement of the average cross section and of its variance is due to the perturbations induced by the electrons. This gives a correlation between the maximum amplitudes of the inter-nuclear oscillational motion and the enhancement factor. It suggests that the chaotic behavior of the electronic motion affects the magnitude of the enhancement factor.PACS numbers: 25.45.-z, 34.10.+x The knowledge of the bare nuclear reaction rates at low energies is essential not only for the understanding of various astrophysical nuclear problems, but also for assessing the effects of the host materials in low energy nuclear fusion reactions in matter. This is currently a subject of great interest in nuclear physics, since Muenster group has reported that the experimental cross sections of the 3 He(d,p) 4 He and of 2 H( 3 He,p) 4 He reactions with gas targets show an increasing enhancement with decreasing bombarding energy with respect to the values obtained by extrapolating from the data at high energies [1]. Many studies attempted to attribute the enhancement of the reaction rate to the screening effects by bound target electrons. In this context one often estimates the screening potential as a constant decrease of the barrier height in the tunneling region through a fit to the data. A puzzle has been that the screening potential obtained by this procedure exceeds the value of the so called adiabatic limit, which is given by the difference of the binding energies of the united atoms and of the target atom and it is theoretically thought to provide the maximum screening potential [2]. Over these several years, the redetermination of the bare cross sections has been proposed theoretically [3] and experimentally [4], using the Trojan Horse Method [5]. The comparison between newly obtained bare cross sections, i.e., astrophysical S-factors, and the cross sections by the direct measurements gives a variety of values for the screening potential. There are already some theoretical studies performed using the time-dependent Hartree-Fock(TDHF) scheme [6,7].In this letter we examine the subject within the constrained molecular dynamics (CoMD) model [8], even in the very low incident energy region not reached experimentally yet. At such very low energies fluctuations are anticipated to play a substantial role. Such fluctuations are beyond the TDHF scheme. Not only TDHF calculations are, by construction, cylindrically symmetric around the beam axis. Such a limitation is not necessarily true in nature and the mean field dynamics could be not correct especially in presence of large fluctuations. Molecular dynamics contains all possible correlations and fluctuations due to the initial conditions(events). For the purpose of treating quantum-mechanical systems like target atoms and molecules, we use classical equatio...