low molecular weight PEO is commonly called as Poly(Ethylene Glycol) (PEG) and it is being used as dispersing agents
for drugs. It is also being used to seed helical segments in proteins. Since PEO has the ability to solvate many alkali
salts, it is used as a host matrix in Solid Polymer Electrolytes (SPE). The conduction mechanism of alkali ions in SPE
is not well understood and this is one the reasons why PEO has come in limelight in recent years. It is beleived that
the mechanism of ion transport in SPE is quite different from that observed in glassy or non-polymeric solid
electrolytes. In glassy electrolytes, the ions hop from one site to other within the rigid host frame. On contrary, it has
been recently demonstrated through Nuclear Magentic Resonance (NMR) experiments that flexibility of the PEO,
which is the "frame-work" in SPE, is the dominant driving force for the ions to move. Although plethora of
experimental and theoretical studies have been carried out on PEO:salt complexes, relatively little is known about
the structure and dynamics of PEO in its pristine form (i.e. crystalline PEO).
Crystal Structure of PEO :
X-ray diffraction studies carried out by Tadokoro et al have suggested that the unit cell of crystalline PEO is
monoclinic and has four PEO chains in it. Each chain adopts a distorted (7/2) helical structure in the crystalline state.
If all the C-O and C-C bonds were at the minimum of the trans (T) and gauche (G) wells of the torsional potential energy
surface respectively, PEO would have taken an uniform helical structure with a TTG sequence. In the reported crystal
structure one finds a wide distribution of the C-O torsional angle suggesting that some of the C-O bonds do not "sit" at
the minimum of the torsional potential thus leading to a distorted helix. This distortion has been attributed to the
strong intermolecular interactions present in the crystal. The following figures show the crystal structure of PEO.
Figure : The crystal structure of PEO
(a) View along the chain axis (b) View perpendicular to
chain axis
Note : These images were generated using a software called xmovie. This software was
originally developed by Prof. Preston Moore of Philadelphia. Recently, additional features
were added to this by Sagar Sen of MSG, JNCASR.
We have investigated the conformational changes associated with the order-disorder transition and the melting of
crystalline PEO. We have used the Parrinello-Rahman version of isothermal-isobaric molecular dynamics simulations
over a wide range temperatures and at ambient pressure. In the following discussions, I may refer to these
simulations as NPT-F simulations. NPT-F stands for Fully flexible box of constant N (Number of particles), P (Pressure),
and T (Temperature). We have observed a new solid phase of PEO much below its melting point. This phase has been
characterized to be partially disordered with a loss of long-range, in-plane orientational correlations between
molecules. Across the premelting transition, the PEO chains continue to be oriented along the c-axis of the crystal with
the retention of short-range helical conformational order. Interesting dynamical changes across this transition resulting
in sliding diffusion and helical oscillations of the backbone have been studied. Thus, our calculations provide a
microscopic insight into molecular mechanism of melting of PEO. The anisotropic solid state sliding diffusion of the
polymer backbone and the reorientational motion of methylene groups observed in the premelting solid phase of PEO
are consistent with the NMR results. To know more about the results pertaining to this work and for the technical
details, readers are requested to refer the following papers:
(1) M. Krishnan and S. Balasubramanian, Chemical Physics Letters, 385, 351-356 (2004).
(2) M. Krishnan and S. Balasubramanian, Journal of Chemical Physics (submitted)
Vibrational Spectra of crystalline PEO :
We have also explored the vibrational dynamics of crystalline PEO using Normal Mode Analysis (NMA). The basic aim
of NMA is to represent the complex vibrations exhibited by the system around its equilibrium well of the potential
energy surface as a linear combination of 3N independent normal modes. This can be achieved by expanding the
total potential energy of the system in terms of the atomic displacements away from the equilibrium configuration.
Let U(q10,q20,...,q3N0) be the total potential energy of the system at equilibrium configuration denoted by a set
{q10,q20,...,q3N0}. Here, the qi0's are the nothing but the components of the position vectors of N atoms present in
the system. Let qi represent a new state of qi0 displaced away from equilibrium position by an amount µi. Thus,
qi=qi0 + µi.
Let U(q1,q2,...,q3N) be the potential energy of the system displaced slightly away from the equilibrium configuration.
Then, one can expand U(q1,q2,...,q3N) around the equilibrium value as follows:
Since the total force acting on the system is zero at equilibrium, the second term becomes zero. Thus, U(q1,q2,...,q3N)
can be written in terms of the second derivatives of U with respect to atomic coordinates. One can define a
mass-weighted Hessian matrix with componets defined as follows:
The Eigen vectors and Eigen values of the mass-weighted Hessian matrix have rich informations about atomic
displacements and the vibrational frequencies of the normal modes. I would like to point out here that obtaining the
elements of the Hessian matrix associated with complex macromolecular interaction potentials is a non-trivial process.
It is relatively easy to derive the analytical expressions for the elements of Hessian matrix arising out of bond
stretching, angle bending, and van der Waals interactions. This is not the case with torsional and Coulombic
interactions. Since the electrostatic interactions long-ranged, use of Ewald summation or other such methods are
inevitable for the MD simulations of systems with such interactions. The Hessian componets due to the real and
reciprocal space parts of the Ewald summation technique have to be calculated properly. We have come up with an
efficient procedure of obtaining the analytical Hessian elements arising out of potentials that are generic to
macromolecular interactions. This method is simple and easy to program. After deriving the expressions for the
Hessian elements, I wrote a full-fledged NMA code to study the vibrational spectra of crystalline PEO. This code was
written in a modular way using Fortran-90 and has subroutines to compare the forces and Hessian components
evaluated analytically against their numerical counterparts. In the following figures, I have shown the atomic
displacements of the representative normal modes of crystalline PEO. Although the calculations were performed for a
crystal consisting of 32 PEO chains, I show one of the chains that exhibit significant atomic displacements for clarity.
Figure 2 : Atomic displacements of crystalline PEO is shown at different vibrational frequencies. The red spheres
are carbons while the white spheres are oxygens. Hydrogen atoms are not shown for clarity in all the
modes except the one with frequency 1244 cm-1. Arrows denote the directions of atomic displacements
and their lengths are scaled up for the purpose of visualization. The numbers provided in the figure are
the vibrational frequencies given in cm-1. Symmetries of the vibrational modes are also given.
To know more about the vibrational spectra, low frequency collective modes of PEO and the systematic procedure of
deriving the Hessian elements, I recommend the readers to refer the following paper.
(1) M. Krishnan and S. Balasubramanian, Physical Review B, 064304-1 to 064304-10 (2003).
A common element of several complex molecular systems is the presence of an alkyl group. It forms a crucial
component of a wide variety of systems ranging from membranes, micelles, self assembled monolayers, and even
molten salts. Thus a study of the phase diagram of linear and branched alkanes continues to be interesting, and
relevant to the understanding of complex phenomena.
In order to understand how n-heptane behaves under high pressure, we have performed atomistic molecular
dynamics simulations in Parrinello-Rahman version of the isothermal-isobaric ensemble. This study was motivated
by recent high pressure spectroscopic experiments on n-heptane that have reported a liquid-solid and the alluded
solid-solid transition with increasing pressure. Thus, the present simulation study aims at understanding these
transitions and provide atomistic level picture of these transitions. Starting from the stabilized solid phase at 300K
and 10 kbar, we have investigated the range of these two transitions by a gradual decrease and increase of pressure
respectively. Although the solid-liquid transition has clear signatures such as the formation of gauche defects along
the molecular backbone, the present model does not show any sign for a first order solid-solid transition at high
pressures. However, we have observed a dynamical arrest of the rotation of the methyl end groups and shifts in the
vibrational features at such high pressures. Our simulations also point to extensive changes in the environment
around the end groups of the alkane as a function of pressure.
The atomic probability density maps of a methyl carbon and of a methyl hydrogen around a given methyl group
obtained from our study at 10 kbar and 70 kbar are shown in the following figures.
The atomic probability density map of a methyl carbon around The atomic probability density map of methyl
a given methyl group is shown up to its second coordination hydrogen around the carbon atom are shown
shell at 10 and 70 kbar. The black lobes denote the probability at 10 and 70 kbar. Probability density values
density of finding neighboring methyl carbons that belong to greater than 0.0002 hydrogen atoms/A3 are
same layer as the central methyl (grey sphere of the molecule shown.
shown), while the grey lobes denote that of adjacent layer
methyl carbons. The rectangular box is the boundary around
the closest eight neighbors. Probability density values greater
than 0.00008 methyl carbons/A3 are shown.
We have also examined the change in vibrational features of this system as a function of pressure. Some of the
vibrational modes obtained from NMA at 10 kbar are shown in the following figure.
Figure : The atomic displacement vectors of a molecule exhibiting large amplitude motion are shown at different
frequencies: 237 cm-1 (methyl rotation); 817 cm-1 (methyl rocking); 1351 cm-1 and 1353 cm-1 (umbrella
motion of methyl group. The pink sphere represents methyl carbon while the green sphere represents
methylene carbon.
Further details can be obtained from our recent paper :
(1) M. Krishnan and S. Balasubramanian, Journal of Physical Chemistry B, 109, 1936-1946 (2005)
The cover art associated with this paper is shown below.
Ultrathin Crystalline films of n-alkanes adsorbed on graphite :
layers at the solid/fluid interface. The packing and orientation of these adsorbed species decide the properties of the
solid layers. The presence of competing interactions, such as adsorbate-adsorbate and substrate-adsorbate
interactions gives rise to a rich phase diagram. The preferential adsorbtion of n-alkanes on graphite is due to the
reason that one of the lattice constants of graphite (2.46 Å) is quite close to the 1-3 distance (the distance between
methylene groups separated by one CH2 group is about 2.54 Å) in a linear alkane. The study of the complex phase
behaviour exhibited by these nearly two dimensional systems has applications to phenomena such as lubrication,
wetting, and spreading.
In what follows, I describe our simulation works on the phase behaviour of such thin molecular films of alkanes
adsorbed on graphite.
Coverage dependent structural changes in n-hexane films on graphite :
We have performed all-atom molecular dynamics simulations of n-hexane on the basal plane of graphite at monolayer
and multilayer coverages. In keeping with experimental data, we find the presence of ordered adsorbed layers both at
monolayer coverage and when the adsorbed layer coexists with excess liquid adsorbate. Using a simulation method
that does not impose any particular periodicity on the adsorbed layer, we quantitatively compare our results to the
results of neutron diffraction experiments and find a structural transition from a uniaxially incommensurate lattice to a
fully commensurate structure on increasing the coverage from a monolayer to a multilayer. The molecular mechanism
associated with this transition has also been studied. The results of this study are summarized in the following figures.
Top view of the time averaged configuration Neutron diffraction intensity as a function of wavevector.
of the n-hexane (a) momolayer and (b) trilayer (a) experimental data (b) simulated trilayer and (c) simulated
on graphite. Molecules in black have at least monolayer. Miller indices corresponding to the reflections are
one intermolecular distance that is less than shown in brackets.
5.0 A. Hydrogens are not shown for clarity.
The complete details of simulation methodology and more results can be obtained from the following papers.
(1) M. Krishnan, S. Balasubramanian, and S. Clarke, Journal of Chemical Physics 118, 5082-5086 (2003).
(2) M. Krishnan, S. Balasubramanian, and S. Clarke, Proceedings of Indian Academy of Sciences (Chemical Sciences)
115, 663-677 (2003).
Melting of ultrathin n-heptane films on graphite :
We have employed constant temperature molecular dynamics simulations to understand the melting behavior of
n-heptane bilayer films physisorbed on the basal plane of graphite. The calculations were performed for three different
systems: an uniaxially commensurate monolayer (UCM), an uniaxially commensurate bilayer (UCB) and a fully
commensurate bilayer (FCB). For FCB, the simulations were carried out at 35 (42 for UCM and 36 for UCB) different
temperatures with an interval of 5K. For bilayer systems, we have observed a layer-by-layer melting of the molecular
film. The calculated intermolecular energies show the FCB structure to be more stable than the UCB structure thus
supporting the conclusions of x-ray and neutron scattering experiments. In FCB, the intermolecular and the substrate
energy of the bottom layer exhibit discontinuities in slope at 290 K, while such an effect was observed at 270 K and
250 K for the bottom layer of UCB and UCM, respectively.The enhancement in the temperature of melting of the bottom
layer relative to bulk melting is in accordance with calorimetric and di raction experiments. The enhancement in the
temperature of melting of the first adsorbed layer in the bilayer system has two contributions. (i) the increased in-plane
density and (ii) the presence of an additional layer. We have explored their roles by performing simulations of a model
system, that of an uniaxially commensurate bilayer. Our results suggest that the enhancement in the melting
temperature of the bilayer system has significant contributions from both these effects. The transition temperatures
observed in our simulations are higher than what has been observed in experiments. This is likely due to the
superheating of the crystalline phase due to the inevitably large heating rates (around 1011 degree/s) employed in
simulations. The difference between the transition temperatures obtained from experiments and simulations could
also be due to the constant volume conditions employed in our work. The use of periodic boundary conditions
necessitates that the simulation box contain an integral number of graphite unit cells, as the external potential needs
to be periodic. This demands that the simulations be carried out in the canonical ensemble (NVT). System size effects
too could alter the transition temperatures. For FCB, in the temperature range between 220K and 290K (190 K and 270 K for UCB and 165 K to 250 K for UCM), we find a fraction of molecules in the bottom layer to undergo rotation about their long-axis, preserving
their all-trans character, in agreement with STM experiments. These molecular rotations precede the melting of the film.
We have identified those rotated molecules and have determined their spatial correlations within the layer, in the
temperature range of 220K and 290K (190 K and 270 K for FCB and 165 K to 250 K for UCM). The distribution of
such rotated molecules in the bottom layer changes systematically from solid-like to liquid-like until melting when all
correlations are lost.
For more details, refer the following paper:
(1) M. Krishnan and S. Balasubramanian, Physical Chemistry Chemical Physics (Submitted)
Fluctuating Charge model of water :
Hydrogen bonding in small water clusters is studied through computer simulation methods using a sophisticated,
empirical model of interaction developed by Rick et al (S W Rick, S J Stuart and B J Berne 1994 J. Chem. Phys. 101 6141)
and others. The model allows for the charges on the interacting sites to fluctuate as a function of time, depending on
their local environment. The charge flow is driven by the difference in the electronegativity of the atoms within the water
molecule, thus effectively mimicking the effects of polarization of the charge density. The potential model is thus
transferable across all phases of water. Using this model, we have obtained the minimum energy structures of water
clusters up to a size of ten. The cluster structures agree well with experimental data. In addition, we are able to
distinctly identify the hydrogens that form hydrogen bonds based on their charges alone, a feature that is not possible
in simulations using fixed charge models. We have also studied the structure of liquid water at ambient conditions
using this fluctuating charge model.
The minimum energy structures of water clusters obtained from these studies are shown below.
For further details, please refer the following article:
(1) M. Krishnan, A. Verma, and S. Balasubramanian, Proceedings of the Indian Academy of Sciences (Chemical Sciences)
113 , 579-590 (2001)