Multipulse dynamic NMR has been employed to study molecular order and dynamics of deuteron (2H) labeled phospholipid membranes. Variation of pulse sequence and pulse separation provides the large number of independent experiments necessary for a proper molecular characterization of the systems. Analysis of these experiments is achieved by employing a density matrix formalism, based on the stochastic Liouville equation. Arbitrary relaxation rates and line shapes of single and multiple quantum transitions are considered. The various 2H NMR experiments of macroscopically unoriented bilayers of 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC), specifically deuterated at the 6- and 14-position of the 2-chain, are faithfully reproduced by the model. Computer simulations provide the orientational distributions and conformations of the hydrocarbon chains and the correlation times of the various motions. In the Lα phase the correlation times τR∥ and τR⊥ for chain rotation and chain fluctuation are of the order of 10−8 s, while trans–gauche isomerization occurs significantly faster (τJ∼10−10 s). At the main transition all chain motions slow down abruptly. Further cooling in the Pβ′ phase first continuously decreases the motions. However, 10 K below the pretransition (hysteresis), there is another abrupt slow down of the chain dynamics. In the Lβ′ phase at T=265 K all three motions occur with correlation times of 10−6 to 10−5 s. Because of higher activation energies, however, intermolecular chain motions freeze out first on the time scale of a particular NMR experiment. Thus, at temperatures T<210 K, trans–gauche isomerization becomes the dominant process. Detection of this motion is possible even at T=168 K, where τJ is of the order of 10−4 s. Arrhenius plots of the various correlation times provide the motional activation energies. Values of 9<EJ<14 kJ/mol for trans–gauche isomerization correspond to the local character of this process. As expected, the activation energies for chain rotation (50<ER∥ <69 kJ/mol) and chain fluctuation (53<ER⊥ <79 kJ/mol) are substantially higher. The correlation times for methyl group rotation form a continuous straight line on the Arrhenius plot throughout the three phases studied, yielding an activation energy of EJ(CD3) =9.9 kJ/mol. Molecular order of the chains is discussed in terms of two parameters SZZ and SZ′Z′, characterizing the orientational order of the chains as a whole and the conformational order at a particular segment. In the Lα phase the hydrocarbon chains are partially disordered (0.44<SZZ <0.6) and melted, exhibiting segmental order parameters of SZ′Z′ (C-6)∼0.75 and SZ′Z′ (C-13)∼0.35, respectively. As expected, conformational order decreases from the central unit to the terminal one (order gradient). The Pβ′ phase exhibits two different chain order parameters of SZZ ∼0.6 and SZZ ∼0.9, indicating heterogeneous chain packing. A unique structural interpretation of this result is not yet possible since the microscopic heterogeneity is compatible with most proposed models. In the Lβ′ phase we find SZZ >0.95, SZ′Z′ (C-6)>0.95, and SZ′Z′ (C-13)>0.9, consistent with highly ordered, fully extended hydrocarbon chains.