The parallel molecular dynamics program

for the simulation of liquid crystalline molecules

with complex topologies


J.Ilnytskyi, M.R.Wilson


Dept. of Chemistry, University of Durham,

South Road, Durham, DH1 3LE, United Kingdom


We describe the parallel molecular dynamics program GBMOLDD.  This program uses the domain decomposition algorithm and is targeted at large-scale simulations of molecular systems (particularly polymers and liquid crystals) composed of both spherically-symmetric and nonspherical sites. The nonspherical sites can be described either by a Gay-Berne potential or by soft repulsive spherocylinders.  The molecules can be of arbitrary topology and the intramolecular forces are described via standard force fields. The applications to the computer simulations of the liquid-crystalline systems with complex topologies are discussed.



Investigations on the interplay of colloids and

nematic liquid crystals under confinement


Evelina B. Kim, Roland Faller, Sylvain Grollau, and Juan J. de Pablo


Department of Chemical Engineering, University of Wisconsin-Madison, Madison, WI


It has been shown experimentally that liquid crystals can be used to amplify surface events such as the specific adsorption of proteins on the substrate of a liquid crystal cell.  These experiments constitute the basis for the development of optical sensors that rely on liquid crystal ordering. 

To study such systems by simulation, we compare the effectiveness of three different techniques for the calculation of chemical potentials. These are the customary test particle method proposed by Widom, the recently proposed density of states Monte Carlo, and the expanded canonical ensemble methods. We apply these in a system of bulk liquid crystal modeled as Gay-Berne mesogens.  The results derived from the three methods agree. The test particle method, however, fails in the nematic regime.

Further, we perform Monte Carlo simulations to study the effects of suspended proteins on the structure of a nematic liquid crystal at a molecular level.  Liquid crystals modeled as repulsive Gay-Berne mesogens interact with the proteins represented as soft repulsive spheres.  We consider the case where all interfaces exhibit homeotropic anchoring conditions.  For small proteins (~5nm), we observe that the particles are surrounded by a Saturn ring disclination line.  We investigate the effect of the proteins size on the defect structure.  Defect structure from the simulations serves as a starting point for continuum elastic theory.  We compare our results with theoretical calculations performed for this system.

Finally, we consider a system where a spherical particle is suspended in a confined nematic liquid crystal. We calculate the potential of mean force between the sphere and a substrate mediated by the liquid crystal by means of Monte Carlo simulations. Three methods are used: canonical Monte Carlo, umbrella sampling, and a novel technique that combines canonical expanded ensemble simulations with the recently proposed density of states formalism. The latter method offers clear advantages over the other ones in that it ensures good sampling of phase space without prior knowledge of the energy landscape of the system.

The resulting potential of mean force, computed as a function of the normal distance between a sphere and a surface, suggests that a sphere is attracted to the surface, even in the absence of attractive molecular interactions.

As a side note, we investigate the effect of the interaction potential on the defect structure.  We perform Monte Carlo simulations of a colloid in a nematic environment under conditions of confinement. The liquid crystal is modeled in two different ways. In the attractive case, LC molecules interact with themselves, the sphere, and the surfaces via a generalized Gay-Berne (GB) potential.  The LCs anchor homeotropicaly to all interfaces.  The surface potential is only dependent on the distance from the surface and the orientation with respect to the plane.  In separate simulations we employ again a purely repulsive potential. Here, the interactions with the surface and the sphere are simplified further to eliminate the dependence on the orientation of the LC particle.  It has been shown that such potentials result in homeotropic anchoring of LC.  It is extremely difficult to find equivalent conditions for such dissimilar systems.  Our criterion for matching the conditions includes density profiles in the bulk-like region of the confined systems.  With a sphere fixed at the center between the surfaces, we found that Saturn Ring defects have very different structures for the two models.  In both cases, a solvation shell is formed close to the sphere.  The order within this shell extends farther away from the colloid for the attractive system.  Also, the order parameter which is close to unity at the surface of the colloid decreases strongly towards the defect core for the attractive system; in the repulsive case this is less pronounced.  A repulsive potential for all interaction is beneficial in several ways: smaller systems may be simulated due to missing cutoff restrictions as well as shorter times may be required to attain equilibrium. Additionally, the results corroborate experiments and theory closer especially in regards to defect structures. The only anchoring that can be reached by the applied repulsive potentials, however, is homeotropic.  In contrast, any type of anchoring is possible with the attractive potential and relative strength of interactions can be fine tuned.



Lattice Boltzman nemato-dynamics

with variable order parameter.


S.V.Lishchuk, C.M.Care, I.Halliday


Materials Research institute, Sheffield Hallam University,

Pond Street, Sheffield, S1 1WB, United Kingdom


A lattice Boltzmann (LB) scheme is described which recovers the equations developed by Qian-Sheng for the hydrodynamics of a nematic liquid crystal with a tensor order parameter. The standard LB scalar density is enhanced to a tensor quantity and the macroscopic momentum, density and tensor order parameter are recovered from appropriate moments of this mesoscopic density. A single lattice Boltzmann equation is used with a direction dependent BGK collision term. Additional forcing terms are used to recover the antisymmetric terms in the stress tensor.  A Chapman Enskog analysis is presented which demonstrates that the Qian-Sheng scheme is recovered provided a lattice with sixth order isotropy is used. The method is validated against analytical results for a number of cases including flow alignment of the order tensor and the Miesowicz viscosities in the presence of an aligning magnetic field. It is also hoped to present results from a LB model of a mixture of a nematic liquid crystal and an isotropic fluid.



Searching for nano-structures of the

liquid crystalline cubic phase


Makoto Yoneya


Japan Science and Technology Corporation,

ERATO Yokoyama Nano-structured Liquid Crystal Project, 

5-9-9 Tokodai, 300-2635, Tsukuba, Japan.


Nano-structures of a thermotropic cubic phase forming liquid crystal molecule, 1,2-bis-[4-n-octyloxy-benzoyl]-hydrazine (BABH8) were studied by molecular dynamics (MD) simulations. The model cubic phase structures were proposed for BABH8 system which are constructed from a locally orientationally ordered “bundle” as a building unit and the cell parameter and the space group (Ia3d) from the recent X-ray results. Stability of the model structure was studied by multi-nanosecond MD simulations. Destabilization of the model structure to the columnar phase and stabilization with enhanced hydrogen-bonding were obtained in these simulations.



Multiple Glassy States:

boundaries of the attractive glass.


Christophe Laforge


Department of chemical engineering,  University of Amsterdam,

Nieuwe Achtergracht, 166, 1018 WV, Amsterdam,  The Netherlands


In this poster, we will show how it is possible to probe experimentaly a crossover between the attractive and the repulsive glass in a  monatomic glass-former.



Insertion of a linear homopolymerin a lamellar phase

of surfactants: a numerical study

with a coarse grain model




Fakultät für Physik, Universität Bielefeld

Postfach 100131, D-33501 Bielefeld, Germany


We study  single polymer chains confined between amphiphilic bilayer by numerical simulations of a coarse grained model.  The model contains solvent particles, amphiphilic molecules, and a single linear homopolymer.


First, we study the binary system (solvent + amphiphiles). At high concentration of amphiphilic molecules, the amphiphiles self-assemble into a lamellar phase: a stack of bilayers. The undulations of the bilayer positions agree with  the theoretical prediction for thermally fluctuating fluid membranes with no surface tension.


This  lamellar matrix is then used to confine a linear homopolymer between two bilayers. The polymer is soluble in the solvent, but not in the inner part of the bilayers.



Realistic Monte Carlo simulation of 5CB at a surface.


Abdon Pijpelink


University of Amsterdam, Department of Chemical Engineering,

Nieuwe Achtergracht 166, 1018 WV, Amsterdam, Netherlands


Many technological applications of liquid crystals, in particular liquid crystal displays, depend on orientational effects caused by an alignment layer on a surface. Therefore, the behaviour of liquid crystals at surfaces has been widely studied. Next to the technological importance, understanding the behaviour of liquid crystals at a surface is also of fundament interest.


Experimental techniques like Scanning Tunnelling Microscopy and Molecular Dynamics simulations of mono- or bilayers of adsorbed liquid crystal molecules at graphite or polymer surfaces have provided insight into the organisation of these molecularly thin layers. These studies are only partially relevant for the case of bulk liquid in contact with a surface since the effect of the liquid crystal on top of the layer adjacent to the surface is largely neglected. The complete interfacial region between surface and bulk has also been simulated, using coarse-grained models, treating a complete liquid crystal molecule as a single anisotropic particle. Unfortunately, these simulations lack the atomistic detail that can be of large influence in the interaction between a surface and a liquid crystal.


We present Monte Carlo studies of large systems of 4-n-pentyl-4\'-cyanobiphenyl (5CB) at simple surfaces. The molecules are modelled realistically using a United Atom (UA) approach, in which the molecular structure of the real molecule is present. Both bulk 5CB as well as 5CB at a surface have been simulated, looking at the structural and orientational ordering of the liquid crystal.








D.L. Cheung 1,2 , S.J. Clark 2 , and M.R. Wilson1


1 Dept of Chemistry, University of Durham, South Road, Durham, DH1 3LE, UK

2 Dept. of Physics, University of Durham, South Road, Duharm, DH1 3LE, UK


 Equilibrium molecular dynamics simulations have been performed for the nematic phase of the mesogen 4-( trans-4-n-pentylcyclohexyl)benzonitrile (PCH5). The simulation was performed using a  harmonic force field with sites representing hydrogen atoms explicitly included. Long range electrostatic interactions were modelled using an Ewald sum. Simulation data has been analysed to calculate important material properties; the rotational viscosity co-efficent, g1, and the flexoelectric co-efficents, es and  eb . 


g1 has been calculated using several methods: using the angular velocity correlation function of the nematic director n [1], the mean-squared displacement of n  [1]  and statistical mechanical methods based on the rotational diffusion co-efficient [2]. The flexoelectric co-efficients were calculated using a linear response formalism [3]  where es  and eb  are related to the correlation between the polarisation and orientational stress. We find good agreement between experimental values and those calculated from simulation. 


[1] S. Sarman, and D. J. Evans,  J. Chem. Phys. (1990)  99, 9021.

[2] A. V. Zahkarov, A. V. Komolkin, and A. Maliniak,  Phys. Rev. E  (1999) 59, 6802.

[3] V. B. Nemtsov, and M. A. Osipov,  Sov. Phys. Crystallogr. (1986)   31, 125.







I. Cacelli, G. Cinacchi, G. Prampolini, A. Tani.


Dipartimento di Chimica, Università di PISA,

Via Risorgimento 35, 56126 PISA, Italy


When dealing with large molecules, as those forming mesophases, their complexity rules out the possibility of employing  “on-the-fly” methods and makes straightforward calculation at an ab initio level computationally very demanding.

To circumvent the problem we have developed an approach based on decomposing the molecule into relevant fragments, whose interaction potentials can be ab initio calculated and employed to reconstruct the potential energy surface (PES) of a pair of the whole molecules.

This resulting PES can be then fitted on known model potential functions at desired level of realism. Here we have applied the method to study the polyphenyls series modelled as strings of quadrupolar Gay-Berne discs.



Ordering in a fluid of hard spherocylinders:

An entropy-based approach


D. Costa + F. Saija *, and P. V. Giaquinta +


+ Istituto Nazionale per la Fisica delta Materia (INFM) and Dipartimento di Fisica

Università degli Studi di Messina, Contrada Papardo, C.P. 50 - 98166 Messina, Italy

* CNR - Lit. per i Processi Chimico-Fisici, sez. Messina, Via La Farina 237 - 98123 Messina, Italy


Hard spherocylinders (cylinders of length L and diameter D capped with two hemispheres at both ends) provide a suitable model to investigate entropy driven phase transitions in real colloidal fluids, composed of rigid rodlike molecules. The phase behaviour of such system, and in particular the onset and relative stability of different liquid-crystalline mesophases, has been recently characterized over the whole L/D range [1].

We have performed extensive Monte Carlo simulations in the range 3 < LjD < 5, and up to L/ D = 20, in order to investigate the confignrational entropy of the system in the context of the multi-particle correlation expansion originally derived by Nettleton and Green [2]. The onset of nematic and smectic order has been analyzed in the one-phase scheme formerly proposed by Giaquinta and Giunta [3]. This scheme is based on the analysis of a quantity, the so-called residual multi-particle entropy (RMPE), that includes the re-summed contributions of all correlations involving more than two particles [3]. The vanishing of the RMPE signals the structural changes which take place in the system.

We have found that the ordering threshold detected using the zero-RMPE condition systematically correlates with the corresponding transition point, whatever the nature of the higher-density phase coexisting with the isotropic fluid. We have then resolved the translational and angular contributions to the configurational entropy of the system; as expected, it is the orientational contribution that drives the transition while the translational contribution plays only a minor role [4].


[1] P. Bolhuis and D. Frenkel, J. Chem. Phys. 106 666 (1997).

[2] R. E. Nettleton and M. S. Green, J. Chem. Phys. 29 1365 (1958).

[3] P. V. Giaquinta and G. Giunta, Physica A 187 145 (1992).

[4] D. Costa, F. Saija, and P. V. Giaquinta, Chem. Phys. Lett. 282 86 (1998).



Simulations of liquid crystal hydrodynamics.


Enzo Orlandini


Department of Physics, University of Padova,

via Marzolo 8, 35121, Padova,



We present a lattice Boltzmann algorithm for liquid crystal hydrodynamics.  The coupling between the tensor order parameter and the flow is treated consistently allowing investigation of a wide range of non-Newtonian flow behavior.  We present results for the effect of hydrodynamics on defect coalescence; of the appearance of the log-rolling and kayaking states in Poiseuille flow; and for banding and coexistence of isotropic and nematic phases under shear. Lattice Boltzmann simulations are used to explore the behavior of liquid crystals subject to Poiseuille flow. In the nematic regime at low shear rates we find two possible steady state configurations of the director field.  The selected state depends on both the shear rate and the history of the sample.  For both director configurations there is clear evidence of shear-thinning, a decrease in the viscosity with increasing shear rate. Moreover, at very high shear rates or when the order parameter is large, the system transforms to a log-rolling state with boundary layers that may exhibit oscillatory behavior. We study the kinetics of the nematic-isotropic transition in a two-dimensional liquid crystal by using a  lattice Boltzmann scheme that couples the tensor order parameter and the flow consistently. The time dependences of the correlation function, energy density, and the number of topological defects are shown to obey dynamic scaling laws with growth exponents that, within the numerical uncertainties, agree with the value 1/2 expected from simple dimensional analysis. We find that these values are not altered by the hydrodynamic flow.



Bilayered Mesophases From Tapered Molecules.


Roberto Berardi, Matteo Ricci, and Claudio Zannoni


Dipartimento di Chimica Fisica e Inorganica, Università di Bologna,

Viale Risorgimento 4, 40136 Bologna, Italy.


The understading of molecular factors that could lead to a specific mesophase organization, is a goal of great interest.  Here we develop a simple molecular model [1] for non--centrosymmetric molecules. With suitable parametrisations of both  shape and interaction part  of this model potential we were able to obtain an interesting mesophase behaviour, i.e. a fluid bilayer structure.

We present results of Monte carlo (MC) computer simulation in the isobaric--isothermal ensemble (NPT) for single component systems of 8196 tapered molecules with three interaction--type parameters.  We show how the combined effect of shape and interaction anisotropy could lead to the desired collective behaviour, and how the interaction anisotropy can play a central role in both the long and short range.



Molecular Dynamics simulations

of spherocylinders interacting via induced electric dipoles


M.Rotunno (1,2), T.Bellini (1), M.Glaser (2), Y.Lansac (2) , F.Mantegazza (3)


(1) Dipartimento di Chimica e Biochimica Medica and INFM,

Università di Milano, Milano, Italy

(2) Department of Physics and Ferroelectric Material Research Center,

University of Colorado,  Boulder Colorado, USA

(3) Dipartimento di Medicina Sperimentale and INFM,

Università di Milano Bicocca, Italy


We present a computer simulation study of the effects of dipolar interactions on phase formation and ordering of molecular  fluids.

The dipole-dipole coupling is often neglected among the driving forces governing the stability and phase transitions of  Nevertheless, in some speci_c cases it has been found that dipole interactions promote the large scale organization of condensed matter phases. In particular, several properties of liquid crystals depend on the molecular dipole strength. While recent computer simulation have studied the effect of permanent dipoles on the molecular order, much less attention has been given to the case of dipoles induced by external fields.

We focus on a system of soft spherocylinders [1] characterized by anisotropic polarizability, which, under the effect of an external electric field, interact throught induced dipole-induced dipole (id-id) coupling. Our model is suitable for molecular dynamics.

The id-id interactions are computed by mean of an iterative loop that determines the dipole’s strength of every molecule as deriving both from the external electric field and from the field generated by the dipoles carried by the other molecules.

We performed NVT simulations for a system of 900 spherocylinders. In the case of aspect ratio L/D = 5, the well known [2] isotropic-nematic transition for uncharged hard spherocylinders is reproduced for soft spherocylinders. For large enought field the transition is replaced by a smooth concentration-dependent nematic order. We find that for L/D = 5 such a dependence is weakly affected by the id-id interactions. A di_erent behavior is found for L/D = 2. In this case the crystal-isotropic transition is replaced by a continuous transition which depends not only on the strength of the coupling of each dipole with the external field, but also on the strength of the interactions between the induced dipoles.

In particular we find the id-id interactions decrease the overall alignment of the particles with the external field.


1. M.P.Allen, Phys. Rev. E, 62, 6706 (2000)

2. S.C.McGrother et al:, J. Chem. Phys. 104, 6755 (1996)






Lorna M. Stimson, Jaroslav M. Ilnytskyi and Mark R. Wilson.


Department of Chemistry, University of Durham, South Road, Durham DH1 3LE, U. K.


This paper describes molecular dynamics and Monte Carlo simulation studies of a flexible liquid crystal polymer and a dendritic liquid crystal system.

The liquid crystal polymer study involves molecular dynamics simulations in the NpT ensemble of a siloxane side-chain liquid crystalline polymer. A hybrid model is employed which combines spherical Lennard-Jones sites, to represent the siloxane polymer backbone and flexible alkyl spacers in the side-chain, and anisotropic Gay-Berne sites to describe the mesogenic moieties. 64 separate molecules are used to represent a polymer melt, which is slowly cooled from the high temperature isotropic phase. Simulations are conducted with and without an external magnetic field present. Below a critical field strength the simulations demonstrate microphase separation of the melt into mesogen rich/backbone rich regions as the system is cooled, with the mesogen rich regions showing orientational ordering of the Gay-Berne particles. Above this critical level of external field these domains line up and we observe a smectic liquid crystal phase. Results are presented for the orientational order, anisotropic structure and dynamics of the LC polymer.

The studies of dendritic liquid crystal system model the first four generations of a carbosilane dendrimer, which have been investigated using a fully atomistic model, employing Monte Carlo techniques. The simulations have been carried out on single molecules using a variable strength mean field potential to represent the influence of surrounding molecules in a nematic phase. The overall molecular shape has been monitored through the moment of inertia spheroid. We find that as a function of mean field strength the molecules change shape from spherical to rod-like.

A simplified model of the same LC dendrimer has also been studied in the presence of a soft repulsive spherocylinder solvent. Here we observe a similar change in the structure of the dendrimer as the solvent is cooled from the isotropic phase into nematic and smectic A phases.



Constant Pressure Molecular Dynamic Simulation Method

for Anisotropic Liquids




JST ERATO Yokoyama Nano-structured Liquid Crystal Project,

TRC 5-9-9 Tokodai, Tsukuba 300-2635 Japan


By introducing an anisotropic factor in the cell dynamics of constant pressure molecular dynamics simulations, we dramatically lessen the artifacts related to cell shapes and overcomes the difficulties of  simulating anisotropic molecules under hydrostatic pressure. The method is especially effective for anisotropic liquids, such as smectic liquid crystals and membranes, however can also be used for crystals and isotropic liquids as well. Crystal-smectic-nematic phase transitions in systems of soft spherocylinders are observed.









JST ERATO Yokoyama Nano-structured Liquid Crystal Project,

TRC 5-9-9 Tokodai, Tsukuba 300-2635 Japan


Some non-amphiphilic thermotropic liquid crystal (LC) molecules have been known as they can form stable Langmuir monolayers on air-water interface. However, formation of such a Langmuir monolayer are not well understood yet except some LC molecules, e.g. cyano-biphenyl derivatives, which have amphiphilic-like characters. We apply a molecular simulation approach to study this subject with special attentions on the LC - water interactions. Molecular dynamics (MD) simulations were done for various terminally alkyl and/or alkoxy substituted azobenzen LC molecules on air-water interface with using the realistic (LC and water) molecular models. The simulation results were compared with those of corresponding amphiphilic modification, i.e. terminal omega-carboxyalkoxy substituted azobenzen. We found that the interaction energetics were not so different between the non-a amphiphilic LC and its amphiphilic modification. Although, increased adhesive (to water layer) energy was confirmed for the amphiphilic modification, cohesive energy between LC molecules were dominant to the adhesion in the both amphiphilic and non-amphiphilic LCs for broad range of area/molecule.









JST ERATO Yokoyama Nano-structured Liquid Crystal Project,

TRC 5-9-9 Tokodai, Tsukuba 300-2635 Japan


Spontaneous optical resolution in a fluid smectic phase, rather than in a crystalline phase, was first observed for a racemic modification of beta-Me-TFMHPOBC. The bent conformation at the double stereogenic part of the (R,S)-beta-Me-TFMHPOBC has been preferentially localized for its fragment molecule by ab initio molecular orbital (MO) calculations and suggested as the origin of its anti ferroelectric property. This strongly fixed bent conformation was pointed out as the important factor for the spontaneous resolution. However, actual relationship between this conformational characteristic and spontaneous resolution has not been clarified. In this study, we try to clarify this relationship with molecular simulation approach. By semi-empirical MO calculations, we found that the molecular bend for the two enantiomers were opposite directions each other strictly depending on their (R,S) and (S,R). It means that the local mirror symmetry at the stereogenic part is magnified to the large molecular bent in this beta-Me-TFMHPOBC. Comparison with the single stereogenic TFMHPOBC and the double stereogenic beta-Me-MHPOBC were also done.