CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. In addition, the CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This paper provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM paper in 1983.

A significant modification to the additive all-atom CHARMM lipid force field (FF) is developed and applied to phospholipid bilayers with both choline and ethanolamine containing head groups and with both saturated and unsaturated aliphatic chains. Motivated by the current CHARMM lipid FF (C27 and C27r) systematically yielding values of the surface area per lipid that are smaller than experimental estimates and gel-like structures of bilayers well above the gel transition temperature, selected torsional, Lennard-Jones and partial atomic charge parameters were modified by targeting both quantum mechanical (QM) and experimental data. QM calculations ranging from high-level ab initio calculations on small molecules to semi-empirical QM studies on a 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayer in combination with experimental thermodynamic data were used as target data for parameter optimization. These changes were tested with simulations of pure bilayers at high hydration of the following six lipids: DPPC, 1,2-dimyristoyl-sn-phosphatidylcholine (DMPC), 1,2-dilauroyl-sn-phosphatidylcholine (DLPC), 1-palmitoyl-2-oleoyl-sn-phosphatidylcholine (POPC), 1,2-dioleoyl-sn-phosphatidylcholine (DOPC), and 1-palmitoyl-2-oleoyl-sn-phosphatidylethanolamine (POPE); simulations of a low hydration DOPC bilayer were also performed. Agreement with experimental surface area is on average within 2%, and the density profiles agree well with neutron and x-ray diffraction experiments. NMR deuterium order parameters (SCD) are well predicted with the new FF, including proper splitting of the SCD for the aliphatic carbon adjacent to the carbonyl for DPPC, POPE, and POPC bilayers. The area compressibility modulus and frequency dependence of 13C NMR relaxation rates of DPPC, and the water distribution of low hydration DOPC bilayers also agree well with experiment. Accordingly, the presented lipid FF, referred to as C36, allows for molecular dynamics simulations to be run in the tensionless ensemble (NPT), and is anticipated to be of utility for simulations of pure lipids systems as well as heterogeneous systems including membrane proteins.

CHARMM-GUI Membrane Builder, http://www.charmm-gui.org/input/membrane, is a web-based user interface designed to interactively build all-atom protein/membrane or membrane-only systems for molecular dynamics simulation through an automated optimized process. In this work, we describe the new features and major improvements in Membrane Builderthat allow users to robustly build realistic biological membrane systems, including (1) addition of new lipid types such as phosphoinositides, cardiolipin, sphingolipids, bacterial lipids, and ergosterol, yielding more than 180 lipid types, (2) enhanced building procedure for lipid packing around protein, (3) reliable algorithm to detect lipid tail penetration to ring structures and protein surface, (4) distance-based algorithm for faster initial ion displacement, (5) CHARMM inputs for P21 image transformation, and (6) NAMD equilibration and production inputs. The robustness of these new features is illustrated by building and simulating a membrane model of the polar and septal regions of E. coli membrane, which contains five lipid types: cardiolipin lipids with two types of acyl chains and phosphatidylethanolamine lipids with three types of acyl chains. It is our hope that CHARMM-GUI Membrane Builder becomes a useful tool for simulation studies to better understand the structure and dynamics of proteins and lipids in realistic biological membrane environments.

We present an all-atom additive empirical force field for the hexopyranose monosaccharide form of glucose and its diastereomers allose, altrose, galactose, gulose, idose, mannose, and talose. The model is developed to be consistent with the CHARMM all-atom biomolecular force fields, and the same parameters are used for all diastereomers, including both the α-and β-anomers of each monosaccharide. The force field is developed in a hierarchical manner and reproduces the gas-phase and condensed-phase properties of small-molecule model compounds corresponding to fragments of pyranose monosaccharides. The resultant parameters are transferred to the full pyranose monosaccharides and additional parameter development is done to achieve a complete hexopyranose monosaccharide force field. Parametrization target data include vibrational frequencies, crystal geometries, solute -water interaction energies, molecular volumes, heats of vaporization, and conformational energies, including those for over 1800 monosaccharide conformations at the MP2/ cc-pVTZ//MP2/6-31G(d) level of theory. Though not targeted during parametrization, free energies of aqueous solvation for the model compounds compare favorably with experimental values. Also well-reproduced are monosaccharide crystal unit cell dimensions and ring pucker, densities of concentrated aqueous glucose systems, and the thermodynamic and dynamic properties of the exocyclic torsion in dilute aqueous systems. The new parameter set expands the CHARMM additive force field to allow for simulation of heterogeneous systems that include hexopyranose monosaccharides in addition to proteins, nucleic acids, and lipids.

We present an extension of the CHARMM hexopyranose monosaccharide additive all-atom force field to enable modeling of glycosidic-linked hexopyranose polysaccharides. The new force field parameters encompass 1→1, 1→2, 1→3, 1→4, and 1→6 hexopyranose glycosidic linkages, as well as O-methylation at the C 1 anomeric carbon, and are developed to be consistent with the CHARMM all-atom biomolecular force fields for proteins, nucleic acids, and lipids. The parameters are developed in a hierarchical fashion using model compounds containing the key atoms in the full carbohydrates, in particular O-methyl-tetrahydropyran and glycosidic-linked dimers consisting of two molecules of tetrahyropyran or one of tetrahydropyran and one of cyclohexane. Target data for parameter optimization include full two-dimensional energy surfaces defined by the Φ/Ψ glycosidic dihedral angles in the disaccharide analogs as determined by quantum mechanical MP2/cc-pVTZ single point energies on MP2/6-31G(d) optimized structures (MP2/cc-pVTZ//MP2/6-31G(d)). In order to achieve balanced, transferable dihedral parameters for the Φ/Ψ glycosidic dihedral angles, surfaces for all possible chiralities at the ring carbon atoms involved in the glycosidic linkages are considered, resulting in over 5000 MP2/cc-pVTZ//MP2/6-31G(d) conformational energies. Also included as target data are vibrational frequencies, pair interaction energies and distances with water molecules, and intramolecular geometries including distortion of the glycosidic valence angle as a function of the glycosidic dihedral angles. The model-compound optimized force field parameters are validated on full disaccharides through comparison of molecular dynamics results to available experimental data. Good agreement is achieved with experiment for a variety of properties including crystal cell parameters and intramolecular geometries, aqueous densities, and aqueous NMR coupling constants associated with the glycosidic linkage. The newly-developed parameters allow for the modeling of linear, branched, and cyclic hexopyranose glycosides both alone and in heterogenous systems including proteins, nucleic acids and/or lipids when combined with existing CHARMM biomolecular force fields.

A revision (C35r) to the CHARMM ether force field is shown to reproduce experimentally observed conformational populations of dimethoxyethane. Molecular dynamics simulations of 9, 18, 27, and 36-mers of polyethylene oxide (PEO) and 27-mers of polyethylene glycol (PEG) in water based on C35r yield a persistence length lambda = 3.7 A, in quantitative agreement with experimentally obtained values of 3.7 A for PEO and 3.8 A for PEG; agreement with experimental values for hydrodynamic radii of comparably sized PEG is also excellent. The exponent upsilon relating the radius of gyration and molecular weight (R(g) proportional, variantM(w)(upsilon)) of PEO from the simulations equals 0.515 +/- 0.023, consistent with experimental observations that low molecular weight PEG behaves as an ideal chain. The shape anisotropy of hydrated PEO is 2.59:1.44:1.00. The dimension of the middle length for each of the polymers nearly equals the hydrodynamic radius R(h)obtained from diffusion measurements in solution. This explains the correspondence of R(h) and R(p), the pore radius of membrane channels: a polymer such as PEG diffuses with its long axis parallel to the membrane channel, and passes through the channel without substantial distortion.

Energies of 119 conformations of normal alkanes from butane to heptane were calculated at approximately the CCSD(T)/cc-pVQZ level. Energies of gauche (g) conformers relative to trans (t) decrease as chain length increases. In what is termed the "positive pentane effect", adjacent gauche conformers of the same sign are stabilized compared to nonadjacent conformers; e.g., for hexane the energies of tgt, tgg, and gtg are 0.600, 0.930, and 1.18 kcal/mol, respectively. Torsional terms in the CHARMM27 (C27) force field were fit to the calculated QM energies to yield a revised potential, C27r. Molecular dynamics simulations of normal alkanes (heptane, decane, tridecane, and pentadecane) with C27r yield higher populations of gauche states, increased transition rates, and improved agreement with experiment as compared to C27. In addition, C27r simulations of a hydrated DPPC lipid bilayer yield improved agreement with the experimental NMR deuterium order parameters for the aliphatic chain ends.

Lipid areas (Aℓ), bilayer area compressibilities (KA), bilayer bending constants (KC), and monolayer spontaneous curvatures (c0) from simulations using the CHARMM36 force field are reported for 12 representative homogenous lipid bilayers. Aℓ (or their surrogate, the average deuterium order parameter in the “plateau region” of the chain) agree very well with experiment, as do the KA. Simulated KC are in near quantitative agreement with vesicle flicker experiments, but are somewhat larger than KC from x-ray, pipette aspiration, and neutron spin echo for saturated lipids. Spontaneous curvatures of bilayer leaflets from the simulations are approximately 30% smaller than experimental values of monolayers in the inverse hexagonal phase.