Volume 140, Issue 10, 14 March 2014

A new bottomup procedure for constructing nonconservative (dissipative and stochastic) interactions for dissipative particle dynamics (DPD) models is described and applied to perform hierarchical coarsegraining of a polar molecular liquid (nitromethane). The distantdependent radial and shear frictions in functionalfree form are derived consistently with a chosen form for conservative interactions by matching twobody forcevelocity and threebody velocityvelocity correlations along the microscopic trajectories of the centroids of Voronoi cells (clusters), which represent the dissipative particles within the DPD description. The Voronoi tessellation is achieved by application of the Kmeans clustering algorithm at regular time intervals. Consistently with a notion of manybody DPD, the conservative interactions are determined through the multiscale coarsegraining (MSCG) method, which naturally implements a pairwise decomposition of the microscopic free energy. A hierarchy of MSCG/DPD models starting with one molecule per Voronoi cell and up to 64 molecules per cell is derived. The radial contribution to the friction appears to be dominant for all models. As the Voronoi cell sizes increase, the dissipative forces rapidly become confined to the first coordination shell. For Voronoi cells of two and more molecules the time dependence of the velocity autocorrelation function becomes monotonic and well reproduced by the respective MSCG/DPD models. A comparative analysis of force and velocity correlations in the atomistic and CG ensembles indicates Markovian behavior with as low as two molecules per dissipative particle. The models with one and two molecules per Voronoi cell yield transport properties (diffusion and shear viscosity) that are in good agreement with the atomistic data. The coarser models produce slower dynamics that can be appreciably attributed to unaccounted dissipation introduced by regular Voronoi repartitioning as well as by larger numerical errors in mapping out the dissipative forces. The framework presented herein can be used to develop computational models of real liquids which are capable of bridging the atomistic and mesoscopic scales.
 COMMUNICATIONS


Communication: Coordination structure of bromide ions associated with hexyltrimethylammonium cations at liquid/liquid interfaces under potentiostatic control as studied by totalreflection Xray absorption fine structure
View Description Hide DescriptionTotalreflection Xray absorption fine structure (TRXAFS) technique was applied for the first time to an interface between two immiscible electrolyte solutions under potentiostatic control. The hydration structure of bromide ions was investigated at polarized 2octanone/water interfaces. TRXAFS spectra at Br Kedge measured in the presence of hexyltrimethylammonium bromide (C6TAB) were slightly modified depending on the Galvani potential difference ( ). The extended Xray absorption fine structure analysis exposed hydration structure changes of bromide ions at the polarized interface. The coordination structure of bromide ions at the interface could be analyzed as compared with bromide ions dissolved in aqueous solution and Br^{−}exchanged resin having quaternary ammonium groups. The results indicated that bromide ions were associated with C6TA^{+} at the polarized interface. The relative contribution of ion association form of bromide ions with quaternary ammonium groups was enhanced at a potential close to the ion transfer of C6TA^{+}, where the interfacial concentration of C6TA^{+} is increased as a function of .
 Top

 ARTICLES

 Theoretical Methods and Algorithms

Effective electron displacements: A tool for timedependent density functional theory computational spectroscopy
View Description Hide DescriptionWe extend our previous definition of the metric Δr for electronic excitations in the framework of the timedependent density functional theory [C. A. Guido, P. Cortona, B. Mennucci, and C. Adamo, J. Chem. Theory Comput.9, 3118 (2013)], by including a measure of the difference of electronic position variances in passing from occupied to virtual orbitals. This new definition, called Γ, permits applications in those situations where the Δrindex is not helpful: transitions in centrosymmetric systems and Rydberg excitations. The Γmetric is then extended by using the Natural Transition Orbitals, thus providing an intuitive picture of how locally the electron density changes during the electronic transitions. Furthermore, the Γ values give insight about the functional performances in reproducing different type of transitions, and allow one to define a “confidence radius” for GGA and hybrid functionals.

Electrostatics of proteins in dielectric solvent continua. I. An accurate and efficient reaction field description
View Description Hide DescriptionWe present a reaction field (RF) method which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of an electrostatics calculation with polarizable molecular mechanics (MM) force fields. The method combines an approach originally suggested by Egwolf and Tavan [J. Chem. Phys.118, 2039 (2003)] with concepts generalizing the Born solution [Z. Phys.1, 45 (1920)] for a solvated ion. First, we derive an exact representation according to which the sources of the RF potential and energy are inducible atomic antipolarization densities and atomic shielding charge distributions. Modeling these atomic densities by Gaussians leads to an approximate representation. Here, the strengths of the Gaussian shielding charge distributions are directly given in terms of the static partial charges as defined, e.g., by standard MM force fields for the various atom types, whereas the strengths of the Gaussian antipolarization densities are calculated by a selfconsistency iteration. The atomic volumes are also described by Gaussians. To account for covalently overlapping atoms, their effective volumes are calculated by another selfconsistency procedure, which guarantees that the dielectric function ε(r) is close to one everywhere inside the protein. The Gaussian widths σ i of the atoms i are parameters of the RF approximation. The remarkable accuracy of the method is demonstrated by comparison with Kirkwood's analytical solution for a spherical protein [J. Chem. Phys.2, 351 (1934)] and with computationally expensive gridbased numerical solutions for simple model systems in dielectric continua including a dipeptide (AcAlaNHMe) as modeled by a standard MM force field. The latter example shows how weakly the RF conformational free energy landscape depends on the parameters σ i . A summarizing discussion highlights the achievements of the new theory and of its approximate solution particularly by comparison with socalled generalized Born methods. A followup paper describes how the method enables Hamiltonian, efficient, and accurate MM molecular dynamics simulations of proteins in dielectric solvent continua.

Electrostatics of proteins in dielectric solvent continua. II. Hamiltonian reaction field dynamics
View Description Hide DescriptionIn Paper I of this work [S. Bauer, G. Mathias, and P. Tavan, J. Chem. Phys.140, 104102 (2014)] we have presented a reaction field (RF) method, which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of polarizable molecular mechanics (MM) force fields. Building upon these results, here we suggest a method for linearly scaling Hamiltonian RF/MM molecular dynamics (MD) simulations, which we call “Hamiltonian dielectric solvent” (HADES). First, we derive analytical expressions for the RF forces acting on the solute atoms. These forces properly account for all those conditions, which have to be selfconsistently fulfilled by RF quantities introduced in Paper I. Next we provide details on the implementation, i.e., we show how our RF approach is combined with a fast multipole method and how the selfconsistency iterations are accelerated by the use of the socalled direct inversion in the iterative subspace. Finally we demonstrate that the method and its implementation enable Hamiltonian, i.e., energy and momentum conserving HADESMD, and compare in a sample application on AcAlaNHMe the HADESMD free energy landscape at 300 K with that obtained in Paper I by scanning of configurations and with one obtained from an explicit solvent simulation.

Multiscale coarsegraining of nonconservative interactions in molecular liquids
View Description Hide DescriptionA new bottomup procedure for constructing nonconservative (dissipative and stochastic) interactions for dissipative particle dynamics (DPD) models is described and applied to perform hierarchical coarsegraining of a polar molecular liquid (nitromethane). The distantdependent radial and shear frictions in functionalfree form are derived consistently with a chosen form for conservative interactions by matching twobody forcevelocity and threebody velocityvelocity correlations along the microscopic trajectories of the centroids of Voronoi cells (clusters), which represent the dissipative particles within the DPD description. The Voronoi tessellation is achieved by application of the Kmeans clustering algorithm at regular time intervals. Consistently with a notion of manybody DPD, the conservative interactions are determined through the multiscale coarsegraining (MSCG) method, which naturally implements a pairwise decomposition of the microscopic free energy. A hierarchy of MSCG/DPD models starting with one molecule per Voronoi cell and up to 64 molecules per cell is derived. The radial contribution to the friction appears to be dominant for all models. As the Voronoi cell sizes increase, the dissipative forces rapidly become confined to the first coordination shell. For Voronoi cells of two and more molecules the time dependence of the velocity autocorrelation function becomes monotonic and well reproduced by the respective MSCG/DPD models. A comparative analysis of force and velocity correlations in the atomistic and CG ensembles indicates Markovian behavior with as low as two molecules per dissipative particle. The models with one and two molecules per Voronoi cell yield transport properties (diffusion and shear viscosity) that are in good agreement with the atomistic data. The coarser models produce slower dynamics that can be appreciably attributed to unaccounted dissipation introduced by regular Voronoi repartitioning as well as by larger numerical errors in mapping out the dissipative forces. The framework presented herein can be used to develop computational models of real liquids which are capable of bridging the atomistic and mesoscopic scales.

Electronic couplings for molecular charge transfer: Benchmarking CDFT, FODFT, and FODFTB against highlevel ab initio calculations
View Description Hide DescriptionWe introduce a database (HAB11) of electronic coupling matrix elements (H ab ) for electron transfer in 11 πconjugated organic homodimer cations. Highlevel ab inito calculations at the multireference configuration interaction MRCI+Q level of theory, nelectron valence state perturbation theory NEVPT2, and (spincomponent scaled) approximate coupled cluster model (SCS)CC2 are reported for this database to assess the performance of three DFT methods of decreasing computational cost, including constrained density functional theory (CDFT), fragmentorbital DFT (FODFT), and selfconsistent charge density functional tightbinding (FODFTB). We find that the CDFT approach in combination with a modified PBE functional containing 50% HartreeFock exchange gives best results for absolute H ab values (mean relative unsigned error = 5.3%) and exponential distance decay constants β (4.3%). CDFT in combination with pure PBE overestimates couplings by 38.7% due to a too diffuse excess charge distribution, whereas the economic FODFT and highly costeffective FODFTB methods underestimate couplings by 37.6% and 42.4%, respectively, due to neglect of interaction between donor and acceptor. The errors are systematic, however, and can be significantly reduced by applying a uniform scaling factor for each method. Applications to dimers outside the database, specifically rotated thiophene dimers and larger acenes up to pentacene, suggests that the same scaling procedure significantly improves the FODFT and FODFTB results for larger πconjugated systems relevant to organic semiconductors and DNA.

Interlayer potential for hexagonal boron nitride
View Description Hide DescriptionA new interlayer forcefield for layered hexagonal boron nitride (hBN) based structures is presented. The forcefield contains three terms representing the interlayer attraction due to dispersive interactions, repulsion due to anisotropic overlaps of electron clouds, and monopolar electrostatic interactions. With appropriate parameterization, the potential is able to simultaneously capture well the binding and lateral sliding energies of planar hBN based dimer systems as well as the interlayer telescoping and rotation of double walled boronnitride nanotubes of different crystallographic orientations. The new potential thus allows for the accurate and efficient modeling and simulation of largescale hBN based layered structures.

Calculation of nuclear spinspin coupling constants using frozen density embedding
View Description Hide DescriptionWe present a method for a subsystembased calculation of indirect nuclear spinspin coupling tensors within the framework of currentspindensityfunctional theory. Our approach is based on the frozendensity embedding scheme within densityfunctional theory and extends a previously reported subsystembased approach for the calculation of nuclear magnetic resonance shielding tensors to magnetic fields which couple not only to orbital but also spin degrees of freedom. This leads to a formulation in which the electron density, the induced paramagnetic current, and the induced spinmagnetization density are calculated separately for the individual subsystems. This is particularly useful for the inclusion of environmental effects in the calculation of nuclear spinspin coupling constants. Neglecting the induced paramagnetic current and spinmagnetization density in the environment due to the magnetic moments of the coupled nuclei leads to a very efficient method in which the computationally expensive response calculation has to be performed only for the subsystem of interest. We show that this approach leads to very good results for the calculation of solventinduced shifts of nuclear spinspin coupling constants in hydrogenbonded systems. Also for systems with stronger interactions, frozendensity embedding performs remarkably well, given the approximate nature of currently available functionals for the nonadditive kinetic energy. As an example we show results for methylmercury halides which exhibit an exceptionally large shift of the onebond coupling constants between ^{199} Hg and ^{13}C upon coordination of dimethylsulfoxide solvent molecules.

Free energy calculations from adaptive molecular dynamics simulations with adiabatic reweighting
View Description Hide DescriptionWe propose an adiabatic reweighting algorithm for computing the free energy along an external parameter from adaptive molecular dynamics simulations. The adaptive bias is estimated using Bayes identity and information from all the sampled configurations. We apply the algorithm to a structural transition in a cluster and to the migration of a crystalline defect along a reaction coordinate. Compared to standard adaptive molecular dynamics, we observe an acceleration of convergence. With the aid of the algorithm, it is also possible to iteratively construct the free energy along the reaction coordinate without having to differentiate the gradient of the reaction coordinate or any biasing potential.

Evaluation of the grandcanonical partition function using expanded WangLandau simulations. III. Impact of combining rules on mixtures properties
View Description Hide DescriptionCombining rules, such as the LorentzBerthelot rules, are routinely used to calculate the thermodynamic properties of mixtures using molecular simulations. Here we extend the expanded WangLandau simulation approach to determine the impact of the combining rules on the value of the partition function of binary systems, and, in turn, on the phase coexistence and thermodynamics of these mixtures. We study various types of mixtures, ranging from systems of rare gases to biologically and technologically relevant mixtures, such as waterurea and watercarbon dioxide. Comparing the simulation results to the experimental data on mixtures of rare gases allows us to rank the performance of combining rules. We find that the widely used LorentzBerthelot rules exhibit the largest deviations from the experimental data, both for the bulk and at coexistence, while the Kong and WaldmanHagler provide much better alternatives. In particular, in the case of aqueous solutions of urea, we show that the use of the LorentzBerthelot rules has a strong impact on the Gibbs free energy of the solute, overshooting the value predicted by the WaldmanHagler rules by 7%. This result emphasizes the importance of the combining rule for the determination of hydration free energies using molecular simulations.

A quasiclassical mapping approach to vibrationally coupled electron transport in molecular junctions
View Description Hide DescriptionWe develop a classical mapping approach suitable to describe vibrationally coupled charge transport in molecular junctions based on the Cartesian mapping for manyelectron systems [B. Li and W. H. Miller, J. Chem. Phys.137, 154107 (2012)]. To properly describe vibrational quantum effects in the transport characteristics, we introduce a simple transformation rewriting the Hamiltonian in terms of occupation numbers and use a binning function to facilitate quantization. The approach provides accurate results for the nonequilibrium Holstein model for a range of bias voltages, vibrational frequencies, and temperatures. It also captures the hallmarks of vibrational quantum effects apparent in steplike structure in the currentvoltage characteristics at low temperatures as well as the phenomenon of FranckCondon blockade.

Calculation of excitation energies from the CC2 linear response theory using Cholesky decomposition
View Description Hide DescriptionA new implementation of the approximate coupled cluster singles and doubles CC2 linear response model is reported. It employs a Cholesky decomposition of the twoelectron integrals that significantly reduces the computational cost and the storage requirements of the method compared to standard implementations. Our algorithm also exploits a partitioning form of the CC2 equations which reduces the dimension of the problem and avoids the storage of doubles amplitudes. We present calculation of excitation energies of benzene using a hierarchy of basis sets and compare the results with conventional CC2 calculations. The reduction of the scaling is evaluated as well as the effect of the Cholesky decomposition parameter on the quality of the results. The new algorithm is used to perform an extrapolation to complete basis set investigation on the spectroscopically interesting benzylallene conformers. A set of calculations on mediumsized molecules is carried out to check the dependence of the accuracy of the results on the decomposition thresholds. Moreover, CC2 singlet excitation energies of the free base porphin are also presented.

Spectroscopic accuracy directly from quantum chemistry: Application to ground and excited states of beryllium dimer
View Description Hide DescriptionWe combine explicit correlation via the canonical transcorrelation approach with the density matrix renormalization group and initiator full configuration interaction quantum Monte Carlo methods to compute a nearexact beryllium dimer curve, without the use of composite methods. In particular, our direct density matrix renormalization group calculations produce a welldepth of D e = 931.2 cm^{−1} which agrees very well with recent experimentally derived estimates D e = 929.7±2 cm^{−1} [J. M. Merritt, V. E. Bondybey, and M. C. Heaven, Science324, 1548 (2009)] and D e = 934.6 cm^{−1} [K. Patkowski, V. Špirko, and K. Szalewicz, Science326, 1382 (2009)], as well the best composite theoretical estimates, D e = 938±15 cm^{−1} [K. Patkowski, R. Podeszwa, and K. Szalewicz, J. Phys. Chem. A111, 12822 (2007)] and D e =935.1±10 cm^{−1} [J. Koput, Phys. Chem. Chem. Phys.13, 20311 (2011)]. Our results suggest possible inaccuracies in the functional form of the potential used at shorter bond lengths to fit the experimental data [J. M. Merritt, V. E. Bondybey, and M. C. Heaven, Science324, 1548 (2009)]. With the density matrix renormalization group we also compute nearexact vertical excitation energies at the equilibrium geometry. These provide nontrivial benchmarks for quantum chemical methods for excited states, and illustrate the surprisingly large error that remains for 1 state with approximate multireference configuration interaction and equationofmotion coupled cluster methods. Overall, we demonstrate that explicitly correlated density matrix renormalization group and initiator full configuration interaction quantum Monte Carlo methods allow us to fully converge to the basis set and correlation limit of the nonrelativistic Schrödinger equation in small molecules.

Exciton dissociation in the presence of phonons: A reduced hierarchy equations of motion approach
View Description Hide DescriptionCombining the reduced hierarchy equations of motion (HEOM) approach with the Wignerfunction formalism, we investigate nonperturbatively exciton dissociation under the influence of a phonon bath in an organic heterojunction. The exciton is modeled by an electronhole pair with the electron moving in the presence of both an external electric field and the Coulomb attraction potential from the hole. In the absence of a phonon bath, calculated HEOM results reproduce those from the OnsagerBraun theory in weak electric fields. In the presence of a phonon bath, substantial deviations from the OnsagerBraun theory are found, signaling phononinduced quantum effects. Furthermore, time evolution of the spatial current distribution is examined, and an initial spike followed by a polarity change of the transient photocurrent have been recovered.

Exploring the topography of the stressmodified energy landscapes of mechanosensitive molecules
View Description Hide DescriptionWe propose a method for computing the activation barrier for chemical reactions involving molecules subjected to mechanical stress. The method avoids reactant and transitionstate saddle optimizations at every force by, instead, solving the differential equations governing the force dependence of the critical points (i.e., minima and saddles) on the system's potential energy surface (PES). As a result, only zeroforce geometry optimization (or, more generally, optimization performed at a single force value) is required by the method. In many cases, minima and transitionstate saddles only exist within a range of forces and disappear beyond a certain critical point. Our method identifies such forceinduced instabilities as points at which one of the Hessian eigenvalues vanishes. We elucidate the nature of those instabilities as fold and cusp catastrophes, where two or three critical points on the forcemodified PES coalesce, and provide a classification of various physically distinct instability scenarios, each illustrated with a concrete chemical example.

Charge asymmetry in the rovibrationally excited HD molecule
View Description Hide DescriptionThe recently developed method for performing allparticle nonBornOppenheimer variational calculations on diatomic molecular systems excited to the first excited rotational state and simultaneously vibrationally excited is employed to study the charge asymmetry and the level lifetimes of the HD molecule. The method uses allparticle explicitly correlated Gaussian functions. The nonlinear parameters of the Gaussians are optimized with the aid of the analytical energy gradient determined with respect to these parameters.
 Advanced Experimental Techniques

Getting a grip on the transverse motion in a Zeeman decelerator
View Description Hide DescriptionZeeman deceleration is an experimental technique in which inhomogeneous, timedependent magnetic fields generated inside an array of solenoid coils are used to manipulate the velocity of a supersonic beam. A 12stage Zeeman decelerator has been built and characterized using hydrogen atoms as a test system. The instrument has several original features including the possibility to replace each deceleration coil individually. In this article, we give a detailed description of the experimental setup, and illustrate its performance. We demonstrate that the overall acceptance in a Zeeman decelerator can be significantly increased with only minor changes to the setup itself. This is achieved by applying a rather low, antiparallel magnetic field in one of the solenoid coils that forms a temporally varying quadrupole field, and improves particle confinement in the transverse direction. The results are reproduced by threedimensional numerical particle trajectory simulations thus allowing for a rigorous analysis of the experimental data. The findings suggest the use of a modified coil configuration to improve transverse focusing during the deceleration process.

Novel twostep laser ablation and ionization mass spectrometry (2SLAIMS) of actorspectator ice layers: Probing chemical composition of D_{2}O ice beneath a H_{2}O ice layer
View Description Hide DescriptionIn this work, we report for the first time successful analysis of organic aromatic analytes imbedded in D2O ices by novel infrared (IR) laser ablation of a layered nonabsorbing D2O ice (spectator) containing the analytes and an ablationactive IRabsorbing H2O ice layer (actor) without the analyte. With these studies we have opened up a new method for the in situ analysis of solids containing analytes when covered with an IR laserabsorbing layer that can be resonantly ablated. This soft ejection method takes advantage of the tenability of twostep infrared laser ablation and ultraviolet laser ionization mass spectrometry, previously demonstrated in this lab to study chemical reactions of polycyclic aromatic hydrocarbons (PAHs) in cryogenic ices. The IR laser pulse tuned to resonantly excite only the upper H2O ice layer (actor) generates a shockwave upon impact. This shockwave penetrates the lower analytecontaining D2O ice layer (spectator, a nonabsorbing ice that cannot be ablated directly with the wavelength of the IR laser employed) and is reflected back, ejecting the contents of the D2O layer into the vacuum where they are intersected by a UV laser for ionization and detection by a timeofflight mass spectrometer. Thus, energy is transmitted from the laserabsorbing actor layer into the nonabsorbing spectator layer resulting its ablation. We found that isotope crosscontamination between layers was negligible. We also did not see any evidence for thermal or collisional chemistry of PAH molecules with H2O molecules in the shockwave. We call this “shockwave mediated surface resonance enhanced subsurface ablation” technique as “twostep laser ablation and ionization mass spectrometry of actorspectator ice layers.” This method has its roots in the wellestablished MALDI (matrix assisted laser desorption and ionization) method. Our method offers more flexibility to optimize both the processes—ablation and ionization. This new technique can thus be potentially employed to undertake in situ analysis of materials imbedded in diverse media, such as cryogenic ices, biological samples, tissues, minerals, etc., by covered with an IRabsorbing laser ablation medium and study the chemical composition and reaction pathways of the analyte in its natural surroundings.
 Atoms, Molecules, and Clusters

Electron correlations and twophoton states in polycyclic aromatic hydrocarbon molecules: A peculiar role of geometry
View Description Hide DescriptionWe present numerical studies of one and twophoton excited states ordering in a number of polycyclic aromatic hydrocarbon molecules: coronene, hexaperihexabenzocoronene, and circumcoronene, all possessing D 6h point group symmetry versus ovalene with D 2h symmetry, within the PariserParrPople model of interacting πelectrons. The calculated energies of the twophoton states as well as their relative twophoton absorption crosssections within the interacting model are qualitatively different from singleparticle descriptions. More remarkably, a peculiar role of molecular geometry is found. The consequence of electron correlations is far stronger for ovalene, where the lowest spinsinglet twophoton state is a quantum superposition of pairs of lowest spin triplet states, as in the linear polyenes. The same is not true for D 6h group hydrocarbons. Our work indicates significant covalent character, in valence bond language, of the ground state, the lowest spin triplet state and a few of the lowest twophoton states in D 2h ovalene but not in those with D 6h symmetry.

D_{3h } [ACE_{3}A]^{−} (E = Al and Ga, A = Si, Ge, Sn, and Pb): A new class of hexatomic monoanionic species with trigonal bipyramidal carbon
View Description Hide DescriptionThe nonclassical trigonal bipyramidal carbon (TBPC) arrangement generally exists as transition states (TSs) in nucleophilic bimolecular substitution (SN2) reactions. Nevertheless, chemists have been curious about whether such a carbon bonding could be stable in equilibrium structures for decades. As the TBPC arrangement was normally realized as cationic species theoretically and experimentally, only one anionic example ([AtC(CN)3At]^{−}) was computationally devised. Herein, we report the design of a new class of anionic TBPC species by using the strategy similar to that for stabilizing the nonclassical planar hypercoordinate carbon. When electron deficient Al and Ga were used as the equatorial ligands, eight D 3h [ACE3A]^{−} (E = Al and Ga, A = Si, Ge, Sn, and Pb) TBPC structures were found to be the energy minima rather than TSs at both the B3LYP and MP2 levels. Remarkably, the energetic results at the CCSD(T) optimization level further identify [GeCAl3Ge]^{−} and [SnCGa3Sn]^{−} even to be the global minima and [SiCAl3Si]^{−} and [GeCGa3Ge]^{−} to be the local minima, only slightly higher than their global minima. The electronic structure analyses reveal that the substantial ionic C–E bonding, the peripheral E–A covalent bonding, and the axial mc2e (multi centertwo electrons) bonding play roles in stabilizing these TBPC structures. The structural simplicity and the high thermodynamic stability suggest that some of these species may be generated and captured in the gas phase. Furthermore, as monoanionic species, their first vertical detachment energies are differentiable from those of their nearest isomers, which would facilitate their characterization via experiments such as the negative ion photoelectron spectroscopy.