Advances in simulations of molecules and materials

CCP5 Annual Meeting 2002

University of Durham, UK

9th - 12th September 2002

ABSTRACTS

INDEX



INVITED TALKS

The prediction of structure with first principles methods

G. Ceder
Department of Materials Science and Engineering, Massachusetts Institute of Technology.

Knowledge of a material's structure is essential to rationalize many of its properties and as such, the prediction of structure is a critical step towards the first principles design of materials. Structural stability at finite temperature is determined by free energy rather than energy, and hence models for the excitations and entropy in the material need to be integrated with quantum mechanics in order to properly describe structural selection.

For crystalline solids a scheme exists whereby reduced Hamiltonians can be defined by successively coarse-graining over faster degrees of freedom. Once this Hamiltonian is parameterized with first principles energy calculations, its free energy can be determined by means of Monte Carlo simulation. It will be shown how this approach has been used successfully to predict the phase diagram of several metallic and oxide systems, in some cases even well before the experimental results had been established.

A radically new approach to make first principles prediction of structure more efficient will also be discussed.

Density functional theory: calculating electronic properties using non-local exchange-correlation functionals

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

Density functional theory is now commonly used to calculate the structural and electronic properties of many materials from first principles. It has been very successful in property prediction over a wide range of applications; lattice parameters, band structures, bulk modulus, molecular systems, liquids to name but a few. However, there are a number of outstanding, but important, properties on which it fails. Many can be traced back to the local nature of the exchange-correlation functional (the only approximation used in the DFT scheme). In this talk I will present a methodology of developing a non-local functional that circumvents many problems of the local functionals and I will also present some electronic structure properties that show its usefulness.

Simulation of polymers - the challenge of wide time and length scales

Julian H.R. Clarke
Chemistry Department, UMIST, Manchester, UK

Dynamic processes in synthetic and biopolymers extend over a very wide range of length and time scales and this presents a particular challenge to computer simulation. Exciting advances in simulation methodology however are now providing information on an increasingly wide range of phenomena. Whilst atomistic level molecular dynamics (MD) simulations gives access to studies of dynamics segmental motion in pure homopolymers the study new 'coarse grain' simulation techniques such as dissipative particle dynamics (DPD) are providing a new approach to simulating polymer solutions, copolymer microphase separation and blend miscibility. In DPD the idea is that the 'particles' of DPD represent groups of real particles (such as colloids or monomers in a polymer) rather than the particles themselves. Hence DPD particles interact through ultrasoft potentials which reflect the finite probability that they will pass through each other! Examples will be discussed of the use of simulations on different time scales to investigate segmental relaxation processes in simple linear polymers, charge transport in ionic heteropolymers, microphase transitions in polymer liquid crystals and long time polymer dynamics in an explicit solvent which show unexpected behaviour in the crossover from the Zimm to Rouse scaling regimes.



Computer simulations of liquid crystal

Claudio Zannoni
Università di Bologna, Dipartimento di Chimica Fisica ed Inorganica and INSTM,
Viale Risorgimento 4, 40136 Bologna, ITALY.

Details of the molecular structure of constituent molecules often have a profound influence on the phase transitions and properties or indeed on the existence of liquid crystals. When a mesogen structure is known, the atomistic approach should allow calculating transition temperatures, molecular organization and phase properties through modelling and computer simulations. This has not much happened in practice until now because of the complexity of mesogenic molecules, although a number of significant atomistic simulations have started to appear [1]. Here we show that reasonable estimates for disordering temperatures now start to become feasible and discuss how this can help our understanding of some classical problems in liquid crystals, such as the odd-even effect [2]. While this progress is encouraging, the atomistic approach requires a well defined input structure, while a problem that has often to be tackled is that of designing molecules that have not yet been synthesized and that can be good candidates to yield mesophases with specific properties of interest (such as ferroelectricity). For this class of problems, the simulation of molecular, rather than atomistic, models allows the identification of some of the key physical features (e.g. aspect ratio and biaxiality, electrostatic moments etc.) that favour a certain collective behaviour, providing hints for the design of novel mesogenic molecules. Here we present a brief summary of the development of computer simulations for systems of particles interacting with Gay-Berne (GB) [3] model potentials and their generalizations [4,5]. We consider in particular the simulation of non chiral ferroelectric liquid crystals from tapered particles [4] and the effect that position and orientation of two outboard dipoles has on the formation of tilted smectic phases [6].

[1] see, e.g. contributions by M. Wilson, M. Glaser, P. Procacci in P. Pasini and C. Zannoni, eds., Advances in the computer simulations of liquid crystals, Kluwer, Dordrecht, 2000.

[2] R. Berardi, L. Muccioli and C. Zannoni, to be published (2002)

[3] J.G. Gay and B.J. Berne, J. Chem. Phys. 74 (1981), 3316.

[4] R. Berardi, M. Ricci and C. Zannoni, Chem. Phys. Chem., 2 (2001).

[5] C. Zannoni,J. Mater. Chem., 11, 2637 (2001)

[6] R. Berardi, S. Orlandi and C. Zannoni, to be published (2002).



Modelling enzyme-catalysed reactions

Adrian J. Mulholland
School of Chemistry, University of Bristol, Bristol BS8 1TS, UK

Computer modelling can provide useful extra insight into the mechanisms of enzyme-catalysed reactions, and complements experimental investigations. In particular, enzyme reactions can be modelled effectively with combined quantum mechanics/molecular mechanics (QM/MM) methods. QM/MM calculations can be carried out with ab initio, or more approximate semiempirical, QM techniques. Applications at both levels will be described. Semiempirical methods can be parameterized to treat a particular reaction accurately, and offer the advantage of being relatively fast, allowing molecular dynamics simulations to be carried out, and free energy profiles to be calculated. We have studied the reactions of a number of enzymes, including flavin-dependent monooxygenases, citrate synthase, and chorismate mutase, with QM/MM techniques. These applications have included tests of the reliability of the QM/MM methods, by comparison with experimental data (e.g. correlation of calculated energy barriers with experimentally observed rate constants), and comparison with higher-level QM calculations. We have also investigated an important model reaction in a glutathione transferase by QM/MM umbrella sampling molecular dynamics simulations, using a specifically parameterized semiempirical model. The calculated free energy profiles, and the ratio of formation of two diastereomeric products, are in good agreement with experiment. Important interactions at the active site have been identified and analysed. Calculation of the effects of mutations suggests factors that may be involved in determining the stereospecificity of different isoenzymes. The QM/MM reaction pathway simulations provide new and detailed insight into the reaction mechanism of this important class of detoxifying enzymes, and illustrate the potential of QM/MM modelling to complement experimental data on enzyme reaction mechanisms.



Beyond the point defect limit: simulation approaches to solid solutions of oxides and highly disordered system

Neil Allan

School of Chemistry, University of Bristol, Cantock's Close, Bristol BS8 1TS

n.l.allan@bris.ac.uk http://dougal.chm.bris.ac.uk/allan

Advances in computational procedures are such that atomistic simulations are now capable of providing detailed information about the structure and thermodynamics of ordered inorganic materials and minerals over a wide range of temperatures and pressures. The vast majority of calculations have studied the properties of perfect crystals, assuming a periodic system in a particular configuration; disorder has largely been investigated theoretically via point defect calculations, which refer only to the dilute limit. Such methods are not readily extended to the solid solutions formed between two similar end-member compounds or disordered systems. Simulations have thus been largely restricted to end-member compounds, excluding many industrially important ceramics and naturally occurring minerals.


This presentation shows explicitly how (i) direct free energy minimisation via lattice quasiharmonic lattice dynamics [1,2] (ii) modified Monte Carlo techniques that allow an efficient sampling of a large number of different configurations can be used for the study of solid solutions, phase diagrams and order-disorder problems in oxides and silicate minerals [2,3]. Explicit account is taken in both of changes in geometry and vibrational effects. We do not resort to the use of a parameterised Hamiltonian that averages out local relaxation and cannot readily be extended to include the effects of lattice vibrations and pressure. Examples include (i) MnO/MgO and CaO/MgO solid solutions. For CaO/MgO we have for the first time calculated the phase diagram for both solid and liquid phases (ii) ordering in a range of spinels and garnets, and pyrope-grossular solid solutions (iii) partitioning of trace elements [4] between oxide melts and solid crystals.



1. M.B. Taylor, G.D. Barrera, N.L. Allan and T.H.K. Barron, Phys. Rev. B56, 14380 (1997); B59, 353, 6742 (1999).

2. N.L. Allan, G.D. Barrera, J.A. Purton, C.E. Sims and M.B. Taylor, Phys. Chem., Chem. Phys. 2, 1099 (2000).

3. M. Yu Lavrentiev, N.L. Allan, G.D. Barrera and J.A. Purton, J. Phys. Chem. 105, 3594 (2001).

4. N.L. Allan, J.D. Blundy, J.A. Purton, M. Yu. Lavrentiev and B.J. Wood, ch. 11 in Solid Solutions in Silicate and Oxide Systems (ed. C.A. Geiger), vol. 3, EMU Notes in Mineralogy, 2001.



CONTRIBUTED TALKS

Ferroelectricity and isotope effects in KDP-type H-bonded crystals

Jorge Kohanoff 1 Sergio Koval, Ricardo Migoni 2 Erio Tosatti 3

  1. Atomistic Simulation Group Queen's University Belfast Belfast BT7 1NN Northern Ireland

  2. Instituto de Fisica Rosario, Rosario, Argentina

  3. SISSA, Trieste, Italy

Based on an accurate first principles description of the energetics in H-bonded KDP, we conduct a first study of nuclear quantum effects and of the changes brought about by deuteration. Correlated tunneling with the pattern of the ferroelectric mode is allowed. It is also found that the main effect of deuteration is a depletion of the proton probability density at the O-H-O bridge center, which in turn weakens its proton-mediated covalent bonding. The ensuing lattice expansion couples selfconsistently with the proton off-centering, thus explaining both the giant isotope effect, and its close connection with geometrical effects.

'Universal' Scaling of the Implantation Depth for 7-Atom Clusters into Graphite

C. F. Sanz-Navarro, R. Smith, S.D. Kenny 1 S. Pratontep, D.J. Kenny, and R. E. Palmer 2

  1. School of Mathematics and Physics, Loughborough University, Loughborough, LE11 3TU, UK.

  2. Nanoscale Physics Research Laboratory, School of Physics and Astronomy, The University of Birmingham, Birmingham, B15 2TT, UK.

Small clusters of atoms deposited on or beneath surfaces have many applications in nanotechnology. It is therefore important to understand the dynamic behaviour of these clusters. The combination of Molecular Dynamics (MD) simulations and experiment can determine how such clusters interact with the surface. Here we report on a "universal" scaling behaviour for the implantation of 7-atom clusters into graphite. It is observed that the implantation depth scales linearly with respect to the cluster momentum according to a similar law for the Si7, Ag7 and Au7 clusters. This implies a constant momentum loss rate per graphite layer almost independent of the velocity and atomic mass of the cluster atoms. Differences between the cluster implantation depth and the actual depth of damage in graphite are understood by analysing the MD simulations. The effect of the experimental oxidative etching technique on the ultimate depth of the pits is also studied by comparison between MD simulations and experiment.

Transport Coefficients of Simple Liquids by Simulation Revisited

D.M. Heyes
Department of Chemistry, University of Surrey, Guildford GU2 7XH.

In this talk I will report the results of Molecular Dynamics simulations of the transport coefficients of hard and soft spheres over a wide density range. The density dependence of the transport coefficients will be compared with the predictions of kinetic theory for hard spheres and semi-empirical expressions. The effects of particle softness on the transport coefficients will be discussed. The relevance of this data to colloidal liquid rheology will be highlighted.

Coil-to-globule transition of isolated homo- dendrimers in solution

Ronan Connolly 1 Edward G. Timoshenko 2 Yuri A. Kuznetsov 3

  1. Theory and Computation Group, Department of Chemistry, University College Dublin, Belfield, Dublin 4, Ireland

  2. Theory and Computation Group, Centre for Synthesis and Chemical Biology, Conway Institute for Biomolecular and Biomedical Research, Department of Chemistry, University College Dublin, Belfield, Dublin 4, Ireland

  3. Centre for High Performance Computing Applications, University College Dublin, Belfield, Dublin 4, Ireland

The regular and well-defined structures of dendritically-branched polymers (dendrimers) provide potential for their application as new advanced materials, superior catalysts, potential drug delivery vehicles amongst other things. In order to better understand the structures of these polymers, conformations of isolated homo-dendrimers of G=1-7 generations with D=1-6 spacers were studied in the good and poor solvents, as well as across the coil-to-globule transition, by means of a version of the Gaussian self-consistent (GSC) method and Monte Carlo (MC) simulations in continuous space based on the same coarse-grained model. This model includes harmonic springs between connected monomers and the pair-wise Lennard-Jones potential with a hard-core repulsion. The scaling law for the dendrimer size, the degrees of bond stretching and steric congestion, as well as the radial density, static structure factor, and asphericity were analysed. It is also confirmed that while smaller dendrimers have a dense core, larger ones develop a hollow domain at intermediate distances from the core.

A parallel molecular dynamics program for the simulation of polymers, dendrimers and other complex molecular topologies.

Jaroslav Ilnytskyi and Mark Wilson
Department of Chemistry, University of Durham, South Road, Durham, DH1 3LE UK

We describe in detail a domain decomposition parallel molecular dynamics program, GBMOLDD. The program is capable of simulating complex molecular topologies composed of both spherical and non-spherical (Gay-Berne and soft spherocylindrical) potentials linked together in arbitrary combinations. NVE, NVT and NPT ensembles are available using a Nose-Hoover thermostat and Hoover barostat with a choice of integrators. The program is targeted at large-scale simulations of polymers, dendrimers and liquid crystalline systems, and shows excellent scalability on a range of parallel machines (including the Cray T3E, shared memory systems and distributed memory clusters). We have used the program to study phase behaviour in model low-molecular weight liquid crystals, liquid crystalline polymers and coarse-grained models of dendritic liquid crystals.

References

[1] A domain decomposition molecular dynamics program for the simulation of flexible molecules with an arbitrary topology of Lennard-Jones and/or Gay-Berne sites

J. M. Ilnytskyi and M. R. Wilson. Comp. Phys. Comm.134 23-32 (2001)

[2] A domain decomposition molecular dynamics program for the simulation of flexible molecules of spherically-symmetrical and nonspherical sites. II. Extension to NVT and NPT ensembles.

J. M. Ilnytskyi and M. R. Wilson. Comp. Phys. Comm. 148 43-58 (2002)



Simulation of Cation Disorder

David J. Willock
Department of Chemistry, Cardiff University, P.O.Box 912, Park Place, Cardiff, CF10 3TB
email: willockdj@Cardiff.ac.uk

Many oxide materials contain a considerable degree of disorder in which crystallographic sites have partial occupancy. The degree to which this disorder will affect the structural and chemical properties of the material can be assessed if the most thermodynamically important structures can be identified.

We have developed an approach to this problem in which a large number of possible structures are simulated and the results stored in a simple encoded form. The software uses the GULP code to calculate the energy of each structure so that the relative energy of each configuration is also retained. As an example we will consider a spinel model of the γ-alumina structure. In this case an exhaustive list of configurations was compiled for a structure consisting of 3 cubic cells with complete cation occupancy of octahedral sites and partial occupancy of tetrahedral. After initial single point energy calculations the results were ranked and the most stable structures optimized in order. Convergence of the calculation was checked by comparing the thermodynamic probability distribution of lattice energy for each set of configurations. The addition of extra structures after 78 137 configurations had been relaxed was found to have little affect on the distribution observed. This calculation gives us direct access to the partition function for the tetrahedral site distributions with long range forces taken into account, allowing the calculation of thermodynamic probabilities. It is shown that good agreement with the structural data on γ-alumina can be obtained from thermodynamic averages of this data.

Applications of this methodology to problems in cation ordering in zeolites will also be discussed.

Recent results in the simulation of liquid crystals

Mike Allen 1 Muataz Al-Barwani 2

  1. University of Warwick

  2. Sultan Qaboos University, Oman

This talk will outline some recent results in the simulation of liquid crystals using variants of the Gay-Berne model. Molecular dynamics simulations using massively-parallel supercomputers and cluster computers have been carried out. The structure of the director field in a nematic liquid crystal around spherical and non-spherical colloid particles has been investigated, and the elastic-mediated interactions between colloidal particles in a nematic solvent determined.



Chiral Indices for Liquid Crystal Molecules

M. P. Neal and M. Solymosi
School of Mathematical and Information Sciences, Coventry University, Priory Street, Coventry, CV1 5FB, U.K

Chirality of liquid crystal molecules has become an important research topic. Recently a range of different approaches have been undertaken to create a better understanding of the concept of chirality in ferroelectric liquid crystals 1,2. Goodby1 examined the dependency of ferroelectricity in liquid crystals as a function of the rotational freedom of the chiral centre and the overall molecular polarity in a homologous series of compounds. A new simulation technique has been published2 that allows an accurate calculation of the helical twisting power (HTP) values for a range of chiral dopant molecules. Those calculated HTP results are relatively consistent with experimental determined values. A newly developed scaling method3 of intrinsic chirality index 4,5 allowed the calculation a chiral measure for a selection of molecular structures 1,2. The normalization inherent in the scaling allows meaningful comparison between molecules with different numbers of atoms. The scaled index is rapid and simple to calculate, the order of seconds on a p.c. It has been successful in predicting the sign and magnitude of the helical twisting power for a range of molecules. We discuss the results for rigid and flexible molecules and relate them to properties of the minimised molecular geometries.

1. J.W. Goodby, J. Mater. Chem., 1991,1, 3, 307-318.

2. M.J. Cook and M.R. Wilson, J. Chem. Phys., 2000,112, 3, 1560-1564.

3. M. Solymosi, R.J. Low, M. Grayson, M.P. Neal, J. Chem. Phys., 2002in press.

4. M.A. Osipov, B.T. Pickup and D.A. Dunmur, Mol. Phys., 1995, 84, 6, 1193-1206.

Towards modelling flexoelectric liquid crystal devices

F.Barmes, C.M.Care, 1 C. Zannoni 2 D.J.Cleaver 1

  1. Materials Research Institute, Sheffield Hallam University

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

The main effects of confinement in liquid crystalline systems, namely surface-induced layering and surface-induced ordering, have been observed experimentally for about two decades [1]. However, molecular-level understanding of these processes is far from complete and the theoretical picture remains largely empirical despite the importance of surface interactions to liquid crystal display technology [2]. Moreover, it seems that the next generation of display will provide an important breakthrough in terms of power consumption by using surface induced switching and flexoelectric mesogens [3,4].
Here, we report results from a series of simulations on confined and unconfined hard particle systems aimed at giving clearer insight into these important processes. For ellipsoidal particles confined in slab geometry we find that bistable anchoring can be generated from a very simple particle-wall interaction. We also show that (tunable) tilted surface anchoring can be generated in a system with purely steric interactions - a result which challenges previous theoretical and simulation work. In the second part of the talk we describe results obtained in the modelling of purely repulsive pear shaped particles [5]. Here we show that such systems have some unexpected phase behaviour, especially when compared with previous simulations of equivalent systems with attractive interactions.



[1] Jerome B. Physics reports 54 391-451 (1991)

[2] Shanks I.A. Contemporary Physics 23 65-91 (1982)

[3] Brown C.V., Towler M.J., Hui V.C and Bryan-Bown G.P. Liquid Crystals 27 233-242 (2000)

[4] Davidson A.J. Mottram N.J. Physical Review E 65 051710 (2002)

[5] Berardi R., Ricci M., Zannoni C. ChemPhysChem 2(7) 443-447 (2001)

Simulations of protein folding and spectroscopy

Jonathan D. Hirst
School of Chemistry University of Nottingham University Park Nottingham NG7 2RD

Recently, we have made significant improvements in the accuracy of calculations of the circular dichroism of proteins from first principles [1]. The quality of these calculations (especially at 220 nm, a key wavelength, where the intensity of the band correlates well with the helical content of polypeptides) has given us confidence to couple such calculations with molecular dynamics simulations of the folding of polypeptides. We use this combined approach to explore the influence of dynamics on the circular dichroism spectroscopy of polypeptides. We apply it to an alanine-based peptide containing an aromatic residue, to probe the effect of the orientation of the aromatic chromophore on the circular dichroism. We analyse a molecular dynamics simulation of the unfolding of myoglobin. Examination of many individual structures allows us to dissect the influence of conformation on the calculated circular dichroism spectra. Our results are aimed at providing a deeper understanding of the optical properties of proteins. An atomic level connection between molecular dynamics simulations and optical spectroscopy is increasingly desirable as theoretical and experimental studies begin to probe protein folding events reliably on the nanosecond timescale.

[1] Besley NA & Hirst JD (1999). J Am Chem Soc, 121, 9636-9644.

Micro-aero-hydrodynamics: new approach to molecular flows in narrow pores and contact line movements on open surfaces

Yu. K. Tovbin
Karpov Institute of Physical Chemistry, 10, ul. Vorontsovo Pole, 105064 Moscow, E-mail: tovbin@cc.nifhi.ac.ru

At present, the dynamics of transport for the mass, momentum and energy in nanoscaled systems for time intervals more than 10-8 -10-7s is usually treated in terms of Navier-Stocks equation. However, this approach is questionable because of the strong fluid inhomogeneity in nanopores as well as contact line movements on open surfaces because of the surface potential, and in case of a menisci formation at the vapor-liquid interface in pores due to molecular interactions. At last time new numerical technique for molecular level investigations on nanoscale temporal and spatial domains have been developed [1-3], which allows to consider various molecular flows in porous systems and for open surfaces during time intervals from 10-11 - 10-9 seconds up to 10-5 seconds. We were using the kinetic theory for the lattice-gas model. The lattice-gas model is the simplest molecular model of the dynamics on the microscopic level. This model takes into account all molecular forces and covers a wide range of fluid densities (from the gaseous up to the liquid state) and temperatures, including the critical range [4,5].

The main advantage of atomic-molecular kinetic equations is that they can be used practically for all times, beginning from microscopic scale (from a characteristic time of atom oscillations) up to macroscopic times, including equilibrium condition of a system [5,6] (kinetic equation for processes in a condensed phase has the Glauber's type equations [7].)

In the lattice-gas model, it is easy to account for correlation effects; therefore this model has a strong advantage over the Boltzmann equation, where only the unary distribution function is used and correlation effects between lateral interacting molecules are absent. The non-equilibrium pair distribution function plays the main role in these new equations, since the form of the equations for the pair distribution function changes at various spatial scales.

The lattice-gas model has some advantages for solving the problem formulated above. First, it provides a self-consistent description of both equilibrium and dynamic properties of the system in a wide range of time scales (starting from the pico-second interval). Second, the calculation based on this model can be performed more rapidly, compared to other numerical methods. Third, the model is in full quantitative agreement with the MD and Monte Carlo simulations in the whole parameter range of the phase diagram, excluding critical points where a qualitative agreement can be achieved, as well as agreement with MD simulation and kinetic theory of inhomogeneous fluids [8]. On basis of the kinetic theory (the Master Equation) the Navier - Stokes type equations were constructed for local equations for the concentration, momentum and energy transfer. Corresponding expressions for concentration dependences of following dissipative factors: selfdiffusion, shear and volume viscosity, as well as thermal conductivity have been obtained. The numerical analysis of the Navier - Stokes type equations are executed by a "relaxation" method of the second order of accuracy on temporary and space coordinates. Using the new approach, we have analyzed the dynamic modes of flows for argon atoms in slit shaped pores in a wide range of densities, including those at which the adsorbed fluid undergoes vapor-liquid separation. At low fluid densities we have obtained a strong anisotropic distribution of fluid on normal to surfaces of pore walls due to adsorption forces. Profiles of speed distributions for flows in pores have been found to differ from the Poiseuille's flow. When the state of the fluid in a nanopore is vapor-like, no viscous flow is observed. When the state of the fluid is liquid-like, the anisotropy of the fluid across the pore was less pronounced, and a volume flow was observed. Using this technique, the dynamic of wetting of a graphite slab by liquid argon was modeled. A detailed study has been performed of the dynamics of the three phase contact line before the steady state regime was established. The velocity of the contact line motion along the Wilhelmy plate is found to decrease as the depth of the solid-fluid potential decreases. In the vicinity of the three-phase contact, the vapor-liquid interface was non-stationary, and formed a foam-type structure. The front of the foam extended along the vapor-liquid and liquid-solid interface. At a lower temperature, the foaming effect was less pronounced.

Acknowledgment. The study was sponsored by the CRDF (project RC2-2214) and RFBR (project 00-03-32153).

References

[1] Yu. K. Tovbin, Rep. Chem. Phys., 21 (2001) 78.

[2] Yu. K. Tovbin, "Molecular fundamentals of microdynamics: molecule transfer in narrow pores", Russ. J. Phys.Chem., 76, (2002) 76.

[3] Yu. K. Tovbin and R.Ya.Tugazakov, Theoretical Foundations of Chemical Engineering, 34 (2000) 99.

[4] T.L.Hill, Statistical Mechanics, McGraw-Hill Book Company, New York, 1956.

[5] Yu. K. Tovbin, Theory of Physical-Chemistry Processes at a Gas-Solid Interface, Mir Publishers & CRC Press Inc., Boca Raton, Florida, 1991.

[6] Yu. K. Tovbin, Progress in Surface Science, 34 (1990) 1.

[7] R.J. Glauber, J.Math. Phys., 4 (1963) 294.

[8] E.Akhmatskaya, B.D.Todd, P.J.Davis, D.J.Evans, K.E.Gubbins, L.A.Pozhar, J. Chem. Phys., 106 (1997) 4684.

Fluid flow in nanopores: Accurate boundary conditions for carbon nanotubes

Vladimir Sokhan, David Nicholson and Nicholas Quirke
Department of Chemistry, Imperial College London, London, UK SW7 2AY

It is well established in gas dynamic studies that the conditions at a surface can influence the reflection of molecules. [1] Two extreme cases were recognised over a century ago, as specular (mirror-like) and diffuse reflection. In the latter, momentum and energy are completely accommodated by the surface, and emission of molecules obeys the cosine law. The whole range of intermediate cases is described by a phenomenological coefficient of Maxwell. [2] Conditions at the wall appear also in the macroscopic hydrodynamic in the form of Dirichlet boundary conditions for the peculiar velocity field. Little could be derived however from these theories about the relationship between the surface properties and reflection conditions. From several factors defining these conditions, the influence of the dynamic state of the solid is probably the less studied. Not surprisingly, most simulation studies of fluid flow in porous solids adopt some measure of ad hoc assumptions in describing fluid solid dynamics. [3]

In this communication, we report the results of simulation of steady-state Poiseuille flow of a simple fluid using realistic many-body potential function for carbon. By analysing temporal profiles of the velocity components of particles colliding with the wall we obtain values of the Maxwell coefficient and for a first time propose boundary conditions for smooth continuum nanotubes such as they equivalent in adsorption, diffusion and fluid flow properties to fully dynamic atomistic model.

  1. E. H. Kennard, Kinetic Theory of Gases (McGraw-Hill, N.Y., 1938).

  2. J. C. Maxwell. Philos. T. Roy. Soc. London 170, 231 (1879).

  3. P. A. Thompson and S. M. Troian. Nature, 389, 360 (1997).

Monte Carlo simulation of microstructure evolution in network forming materials

V. M. Burlakov
Department of Materials, University of Oxford, Parks Road, Oxford, OX1 3PH, UK

Random network model with coordination defects is developed and implemented in simulations of vapour deposition and post annealing of amorphous nonstoichiometric silica [1]. The main feature of the model is that it explicitly takes into account dangling bonds (DB), which play central role in the deposition process. DB migration and annihilation account for the microstructure evolution of the network material during its deposition and post annealing. To cover a bigger time scale of the microstructure evolution a pair bond switching mechanism [2] is used besides the DB related ones to generate new atomic configurations. Using the developed model the porosity of silica and a-Si has been studied. We found that silica possesses an order of magnitude higher porosity due to presence of oxygen atoms, which act as the surfactant lowering the pore surfaces energy. The pore growth kinetics in silica upon simulated annealing has been also studied. We observed continuous pore growth upon annealing, and found the process resembling diffusion driven chemical reaction.

References

[1] V. M. Burlakov, G. A. D. Briggs, A. P. Sutton, Y. Tsukahara, Phys. Rev. Lett. 86, 3052 (2001).

[2] F. Wooten, K. Winer, D. Weaire, Phys. Rev. Lett., 54, 1392 (1985).



Molecular dynamics simulations of surface contact interaction studies

C.W Yong, 1 W. Smith 1 K. Kendall 2

  1. CCLRC, Daresbury Laboratory, Daresbury, Warrington WA4 4AD, UK

  2. School of Chemical Engineering, The University of Birmingham, Edgbaston, Birmingham B15 2TT, UK

The study of two surface contact is important for the understanding of friction. Friction is not a single phenomenon; many processes contribute to the dissipation of bulk kinetic energy. The coefficient of friction, which is introduced to describe contact mechanics in many theoretical treatments and computer simulations, are at best phenomenologically derived and without sound scientific justification at the atomistic level. It is therefore important to investigate systematically the atomistic factors underlying friction and to establish links from these to the macroscopically observed frictional phenomena.

Friction only occurs whenever two bodies come into contact. The aim of this talk is to provide mechanistic view of surface contacts, in atomistic details, of some simple crystalline materials such as MgO and NaCl. As the materials were brought together and subsequently retracted, a rich variety of phenomena has been observed such as 'jump-to-contact', hysteresis, atomic dislocation and neck formation.

Simulation of mineral-melt partition coefficients

Mikhail Yu. Lavrentiev and Neil L. Allan
School of Chemistry, University of Bristol, Cantock's Close, Bristol BS8 1TS

Partition coefficients, defined as the ratio of the concentrations of an impurity element in coexisting liquid and solid host matrices, are very important in geology. We present a new method of calculating partition coefficients using semigrand ensemble Monte Carlo simulations. The results obtained for partitioning behaviour in CaO and other systems agree well with analytical lattice strain models and demonstrate the wide applicability of the method.

Computer Simulations of the Electrolytes and the Interface between Two Immiscible Electrolytes (ITIES) in non-Euclidean Geometries.

Shabnam Hanassab and J. VanderNoot

Chemistry Department, Queen Mary University of London Mile End Road, London E1 4NS, Great Britain.
E-Mails: S.Hanassab@qmul.ac.uk and tjv@Cognitrix.com
Website: www.Qmul.ac.uk/~ca9142/

Canonical ensemble (NVT) Monte Carlo simulations of electrolytes and the interface between two immiscible electrolytes (ITIES) on the 3D surface of a 4D hypersphere will be presented. A hypersphere is a finite, self-contained and isotropic geometry, which is well adapted for the simulations of coulombic systems.

Comparison of our simulations results of basic electrolytes with those of 3D Euclidean simulations reported in the literature showed a very good agreement.

For the first time we present our most recent results of simulations of the ITIES in the hyperspherical geometry. We studied the position of the ions in each phase and concluded that the solvent permittivity of one phase has a direct effect on the distribution of the ions in the neighbouring phase. However, this effect is not the result of the interaction of the ions with the ions in the other phase, but as a result of the interactions of the image forces of the ions in the aqueous and in the organic phase. Our results indicate that as polarity of the organic solvent is lowered the ions in the aqueous phase are repelled from the interface, whereas the ions in the oil phase are attracted to the interface. A consequence of this is that the ITIES cannot be treated as two back to back electric double layers.

It is also possible to apply an external electric field across the ITIES simulation system. This application resulted in separation of the ions in the system. However, we found out that there is higher degree of separation of ions in the organic phase than in the aqueous phase. From the theory of electrostatics, the field within the organic phase (lower permittivity) can even be higher than the externally applied electric field.

In conclusion the 3D non-Euclidean geometry is very efficient and useful in simulations of the electrolytic systems and interfaces.

Graphic1



Solid-fluid boundary from molecular simulation: Calculation of the melting point of NaCl

Jamshed Anwar 1 Daan Frenkel 2 Massimo G. Noro 3

  1. Computational Pharmaceutical Sciences, Department Of Pharmacy, King’s College London, Franklin-Wilkins Building, 150 Stamford Street, London SE1 9NN, UK.

  2. Computational Physics, FOM Institute for Atomic and Molecular Physics Kruislaan 407, 1098 SJ Amsterdam, The Netherlands

  3. Massimo G. Noro Unilever Research Port Sunlight, Quarry Road East, Bebbington, Wirral CH6 33JW, UK.

We report a numerical calculation of the melting point of NaCl. The solid-liquid transition was located by determining the point where the chemical potentials of the solid and liquid phases intersect. To compute these chemical potentials, we made use of free energy calculations. For the solid phase the free energy was determined by thermodynamic integration from the Einstein crystal. For the liquid phase two distinct approaches were employed: one based on particle insertion and growth using the Kirkwood coupling parameter, and the other involving thermodynamic integration of the NaCl liquid to a Lennard Jones fluid. The latter approach was found to be significantly more accurate. The coexistence point at 1074K was characterised by a pressure of - 30 ± 40 MPa and a chemical potential of - 97.9 ± 0.2 kT. This result is remarkably good as the error bounds on the pressure enclose the expected coexistence pressure of about 0.1MPa (ambient). Using the Clausius Clapyron relation, we estimate that dP/dT » 3 MPa/K. This yields a melting point of 1064K±14K at ambient pressure, which encompasses the quoted range for the experimental melting point (1072.45K-1074.4K). The good agreement with the experimental melting-point data provides additional evidence that the Tosi-Fumi model for NaCl is quite accurate. Our study illustrates that the melting point of an ionic system can calculated accurately by employing a judicious combination of free energy techniques. The techniques used in this work can be directly extended to more complex, charged systems.



POSTER ABSTRACTS

Constraint Method for Deriving Non-equilibrium Molecular Dynamics Equations of Motion

Toni Galea and Phil Attard
Ian Wark Research Institute, University of South Australia, Mawson Lakes, SA 5095 Australia

A non-equilibrium molecular dynamics technique is described that is in theory capable of deriving generalised equations of motion for any type of flow. The method is based on constraining the motion of a system in phase space to an appropriate hypersurface. The method is derived with a finite time step specifically for computer implementation, which ensures that the constraints are satisfied at every time step, and hence there is no need for periodic rescaling as with other constraint techniques. This method has the advantage that there is minimal perturbation of the system for two reasons. First the least constraint force is determined in a way similar to Gauss’ principle, and second the constraining factors can be chosen such that the flow can be achieved with the least effect on the system. The method is demonstrated with isokinetic shear flow, in bulk and slit geometries, which illustrates the flexibility of the method. Results for the shear viscosity are in agreement with previously published results.

Ghost Interface Method for the Surface Tension

Michael Moody and Phil Attard
Ian Wark Research Institute, University of South Australia, Mawson Lakes Boulevard, Mawson Lakes, South Australia 5095 Australia

A new formally exact expression for the surface tension has been derived that is suited for computer simulations. This new route to the surface tension is analogous to Widom's expression for the chemical potential which is conceptually different to the traditional Kirkwood-Buff virial route. The method invokes the average interaction between the two well defined uncoupled phases (virtual interface). Besides the planar interface of the liquid-vapour coexisting phase, the technique can be applied to curved interfaces and also to interfaces away from coexistence. Monte Carlo simulation results will be given for the planar interface of a Lennard-Jones fluid.



Molecular Dynamics Simulation of Carbosilane and Poly(propylene imine) Dendrimers: Effect of Explicit Solvent

M.A. Mazo 1 M.Yu. Shamaev 1 E.B. Gusarova 1 N.S. Perov 2 and N.K. Balabaev 3

  1. Institute of Chemical Physics Russian Academy of Science, Kosygin str. 4, Moscow, 119991 Russia
    mazo@polymer.chph.ras.ru;

  2. Institute of Synthetic Polymeric Materials Russian Academy of Science, Profsouznaja str. 70, Moscow, 117393 Russia
    pero@ispm.ru;

  3. Institute of Mathematical Problems of Biology Russian Academy of Science, Pushchino, Moscow Region, 142292, Russia
    balabaev@impb.psn.ru

Series of computer simulation of dendrimers (see [1]) essentially advanced the understanding of spatial structure and molecular mobility of these unusual polymeric structures. However, numerical experiments were generally performed without explicit account of solvent molecules: a solvent influence was simulated by variation of the potential parameters of van der Waals interactions [2,3]. Recently the behaviour of the dendrimer in the explicit solvent has been considered for Meijer dendrimer box [4] and for a crude bead-spring model of dendrimer [5]. In this work, we study the effect of an explicit solvent on the structural and dynamic properties in the framework of a detailed molecular model of two dendrimers of close topology.

Simulations were been carried out for carbosilane and poly(propylene imine) dendrimers of 5th generation under periodic boundary conditions at different temperature. Calculations have been accomplished using the AMBER force field in the united atom approximation. Lennard-Jones particles have been considered as the solvent molecules with potential parameters corresponded to CCl4. There was one molecule of the dendrimer in each calculation cell, and cell's size was large enough to exclude interaction between dendrimers.

The internal organization of dendrimers was analyzed (gyration radii and principal moments of inertia, density distributions for both dendrimer and solvent atoms) as well as dynamic properties of dendrimers (frequencies of conformational transitions, mobility of different center of branching, rotation and translation diffusions), and mobility of solvent molecules. The results obtained were compared with the data of calculations made with implicit account of the solvent and experimental data.

This work was supported by NWO (project no. 99 005 725) and INTAS (project no. 00-00712).

References

  1. "Dendrimers III: Design, Dimension, Function". Topics in Current Chemistry, (Springer, 2001) 212

  2. "Molecular Dynamics Study of Dendrimer Molecules in Solvents of Varying Quality", M. Murat and G.S. Grest, Macromolecules 29, 1278-1285 (1996).

  3. "The influence of the chemical structure of terminal fragments on the spatial-dynamic organization of dendrimers", M.A. Mazo, N.S. Perov, E.B. Gusarova, P.A. Zhilin and N.K. Balabaev, Russ. J. Phys. Chem. 74, Suppl. 1, S52-S58 (2000).

  4. "Effect of repeat unit flexibility on dendrimer conformation as studied atomistic molecular dynamics simulations", C.B. Gorman and J.C. Smith, Polymer 41, 675-683 (2000).

  5. "Statics and dynamics of model dendrimers as studied by molecular dynamics simulations", K. Karatasos, D.B. Adolf and G.R. Davies, J. Chem. Phys. 115, 5310-5318 (2001).

Molecular Dynamics Simulation of Local Fluid Mobility in Micropores

M.A. Mazo 1, A.B. Rabinovich 2, Yu.K. Tovbin 2

  1. Institute of Chemical Physics Russian Academy of Science, Kosygin str. 4, Moscow, 119991 Russia,
    mazo@polymer.chph.ras.ru

  2. Karpov Institute of Physical Chemistry, 10, ul. Vorontsovo Pole, 105064, Moscow;
    tovbin@cc.nifhi.ac.ru

The lattice-gas model (LGM) was successfully applied to the description of mass transfer processes in highly heterogeneous media in last years [1,2]. Considerable attention has been focused on application of this model to the calculation of the basic kinetic coefficients (self-diffusion coefficients, heat conductivity, value and shear viscosity) in narrow pores [3,4]. Comparison of the values of kinetic coefficients obtained by the LGM, by molecular dynamics simulation (MD) and by the kinetic theory of a liquid demonstrated satisfactory agreement [5]. However, fields of this method application as well as accuracy of results obtained are not essentially up to now. In particular, it is no clear, from what time intervals the lattice-gas model is applicable, and how strong is the influence of density and temperature on the model parameters.

A molecular dynamics simulation of Lennard-Jones (LJ) fluid in narrow slit pores and a calculation of same system in the framework of LGM were carried out for the study these problems. The fluid was confined between parallel structured walls [6]. All fluid and wall particles are interacting to each other with cut off LJ (12-6) potential with σ -ε parameters. The masses and dimensions σ of all atoms (wall and fluid) were equal, but fluid-wall and fluid-fluid particle interaction was changed from repulsive to high attractive potential (εwf / ε ff = 9.24). Two pore widths were studied, 5.5 σ and 13 σ , which allows the fluid to form 6 and 12 layers respectively. The fluid temperature was varied over a wide range from freezing to super critical temperatures.

The obtained results demonstrate that the lattice-gas model can be a convenient interpolation tool for quantitative description of dynamics of particles redistribution between layers and the longitudinal self-diffusion coefficients. Comparison of MD simulation results and LGM calculations allows us to determine the parameters of the lattice gas model and to use it for LGM calculations at large time scales.

References

  1. Yu.K.Tovbin, "Molecular fundamentals of microdynamics: molecule transfer in narrow pores", Russ. J. Phys.Chem., 76, (2002) 76-85.

  2. Yu.K. Tovbin, R.Ya. Tugazakov, and N.V. Komarov, "Micro-Hydrodynamic Approach to Molecular Flows in Narrow Pores", Proceeding of 7th International Conference on Fundamentals of Adsorption , Nagasaki, Japan, May 20-25, 2001, 1007-1014.

  3. Yu.K.Tovbin and N.F.Vasyutkin, "Concentration dependences of adsorbate self-diffusion and shear viscosity in narrow slit-like pores", Russian. Chem.Bull., 50 (2001) 1496-1502.

  4. Yu.K. Tovbin, E.E. Gvozdeva, and D.V. Yeremich, "Dynamical characteristics of adsorbate in narrow cylindrical pores", Russ. J. Phys.Chem. (submitted).

  5. E. Akhmatskaya, B.D. Todd, P.J. Davis, D.J. Evans, K.E. Gubbins, and L.A. Pozhar, "A study of viscosity inhomogeneity in porous media", J. Chem. Phys., 106 (1997) 4684-4695.

  6. A.A. Berlin, M.A. Mazo, N.K. Balabaev, Yu.K. Tovbin, "Single-Component and Two-Component Mixtures in Slit-Like Pores: Molecular Dynamics Simulation", Proceeding of 7th International Conference on Fundamentals of Adsorption , Nagasaki, Japan, May 20-25, 2001, 402-409.



Heterogeneous aggregation in binary colloids

K. W. Yu and Andrew C. T. Wong
Department of Physics, Chinese University of Hong Kong, Shatin, NT Hong Kong

Molecular dynamics (MD) simulation has been employed to study the nonequilibrium structure formation of two types of particles in a colloidal suspension, driven by type-dependent forces. We examined the time evolution of structure formation as well as the structural properties of the resulting aggregation by studying the radial distribution function (RDF). The resulting aggregation is well described by a binary colloidal gelation. We compared the structural properties to those for one type of particles. From the MD results, it is evident that there are significant differences between the RDF's of the two cases. Moreover, we found that the average coordination number is generally larger in the monodisperse case for all area fractions considered. Thus, by means of heterogeneous aggregation, it is possible to obtain a wide variety of structures while more close-packed structures are formed for monodisperse colloidal aggregation.

Computer Simulations of Electrolytic Systems in non-Euclidean Geometries.

Shabnam Hanassab and J. VanderNoot

Chemistry Department, Queen Mary University of London Mile End Road, London E1 4NS, Great Britain.
E-Mails: S.Hanassab@qmul.ac.uk and tjv@Cognitrix.com
Website: www.Qmul.ac.uk/~ca9142/

Canonical ensemble (NVT) Monte Carlo simulations of electrolytes and the interface between two immiscible electrolytes (ITIES) on the 3D surface of a 4D hypersphere will be presented. A hypersphere is a finite, self-contained and isotropic geometry, which is well adapted for the simulations of the coulombic systems.

Comparison of our results from simulations of basic electrolytes with those of 3D Euclidean simulations reported in the literature showed a very good agreement, certifying that spherical boundary conditions could be used as an alternative to Ewald summation or reaction field method in Euclidean geometries.

For the first time we present our most recent results for the simulations of the ITIES, accompanied by applying an external field across the system, in non-Euclidean geometries. The effect of the image forces and the polarity of the organic phase on the position and the thermodynamic properties of the system, are studied and explained in detail.

In conclusion we will confirm that hyperspherical geometry can be used effectively in modelling long-range interactions in computer simulations.

Graphic1



Chiral Discrimination by Chemical Force Microscopy: Simulations

Esther Haines, T. Rayment, J. Goodman, and C. Abell
University Chemical Laboratory, Lensfield Road, Cambridge, CB2 1EW

Atomic Force Microscopy (AFM) provides a means of measuring the normal adhesion force and the lateral, or friction, forces between the tip and a sample at length scales less than 100nm. In Chemical Force Microscopy (CFM) the tip and sample are functionalized, for example, by coating them with self-assembled monolayers, so that AFM measures the forces between the functional groups. Thus, CFM provides a way of mapping the spatial distribution of different chemical species on a surface (ref.1). Experiments have shown that this technique can be used to distinguish between enantiomers (ref 2). We have been carrying out simulations to help answer some of the questions that arise, for example, does chiral discrimination require inter-penetration of the surfaces. Both single-molecule simulations and simulations of the interactions between surfaces show differences, depending on the combination of enantiomers.

Monte Carlo simulations of krypton using ab initio potentials

Afshin E. Nasrabad and U. K. Deiters
Physical Chemistry Institute, University of Cologne, Luxemburger Str. 116, 50939 Cologne, Germany

Gibbs ensemble Monte Carlo simulations were performed to calculate the vapor-liquid coexistence properties of pure krypton with a new ab initio potential. The CCSD(T) level of theory and two correlation consistent basis sets were used and an extrapolation method was applied to obtain the basis set limit of the interaction energies. The influence of the three-body potentials in thermodynamic properties of phases were implemented in the Axilrod-Teller triple-dipole form. The simulations results show that the phase behavior of krypton is well described by the ab initio pair potential plus Axilrod-Teller three-body term and the overall agreement of theory with experimental data is good.

Effects of three-body interactions in vapor-liquid phase coexistence of krypton-argon mixture: a Monte Carlo simulation study

Afshin E. Nasrabad and U. K. Deiters
Physical Chemistry Institute, University of Cologne, Luxemburger Str. 116, 50939 Cologne, Germany

This work is in a series of studies which aim at evaluating the adequacy of ab initio pair potentials plus Axilrod-Teller three-body interactions for description of macroscopic properties of pure and binary mixture of noble gases. Counterpoise-corrected ab initio energies resulting from the CCSD(T) level of theory and two correlation consistent basis sets were used to obtain the basis set limit of the interaction energies using an extrapolation method. The resulting potential energy curves were used in the Gibbs ensemble Monte Carlo simulations to calculate the vapor-liquid coexistence for krypton-argon mixture. Calculations are reported in four different isotherms, 193.15, 177.38, 163.15 and 158.15 Kelvins. Comparison of the simulation results and experimental data shows that the inclusion of three-body interactions to the potential leads to a good agreement for all isotherms.



Equilibrium properties of neon by global simulation and an equation of state

A. E. Nasrabad, R. Laghaei and U. K. Deiters
Institut für Physikalische Chemie, Universität zu Köln, Luxemburger Str. 116, 50939

Vapor-liquid equilibria of neon were calculated by Gibbs ensemble Monte Carlo simulations using an ab initio pair potential (global simulation). Simulations were performed from temperatures just above the triple point up to the critical region. Effects of three-body interactions on the phase behavior are reported. Comparison of experimental data with simulation results shows that addition of Axilrod-Teller three-body interaction to the total energy can improve the results, however, there are still large discripancies especially in the liquid branch. It is shown that these discrepancies are mainly due to the quantum effects, and that they can be accounted for by a newly developed correction function.

Modelling template-oriented crystal growth

J H Harding (presenter) and D M Duffy (lead author)
Dept. Physics & Astronomy University College London

The interfaces between the surfaces of calcite crystals and Langmuir monolayers of octadecanoic acid have been modelled using molecular dynamics techniques. All of the seven interfaces modelled were found to be strongly bound, with low interfacial energies and high densities of hydrogen bonds across the interface. The (001) Ca terminated surface was found to have the lowest interfacial energy. The adsorbed monolayers were found to have a significant effect on the surface structure for some crystallographic orientations.

Computation of second virial coefficients of electro-optic phenomena.

V W Couling
School of Chemical and Physical Sciences, University of Natal, Private Bag X01, Pietermaritzburg 3209, South Africa

Theoretical studies of the effects of molecular interactions on the bulk electro-optic properties of fluids have recently been extended to include molecules with symmetry lower than linear. This has been achieved by invoking the powerful algebraic manipulation package Macsyma to perform the tensor analysis and hence to obtain expressions for the virial coefficients. The picture which emerges from the molecular phenomena studied thus far, namely depolarized Rayleigh light-scattering and electric-field-induced birefringence, is that comprehensive molecular-tensor theories of these effects explain the observed phenomena adequately.



Computation and measurement of second Kerr-effect virial coefficients of dimethyl ether.

T J Sono and V W Couling
School of Chemical and Physical Sciences, University of Natal, Private Bag X01, Pietermaritzburg 3209, South Africa

The molecular theory of the second Kerr-effect virial coefficient Bk, describing the effects of interacting pairs of molecules to the electro-optic Kerr effect, is presented for interacting asymmetric-top molecules. Results of the numerical computations of Bk for dimethyl ether are given. The calculated Bk values generally lie within the uncertainty limits of the available measured data over their full range of temperatures. We have used our recently-commissioned apparatus to take precise new measurements of the Kerr effect for dimethyl ether, and the measured data are in good agreement with the findings of our molecular model. It is also possible to deduce anisotropic molecular polarizabilities and molecular hyperpolarizability tensors from the measured Kerr-effect data. The polarizabilities are important properties in many areas of physics and chemistry in that they describe the interaction of a molecular charge distribution with an electric field, and so are responsible for a variety of different effects ranging from refraction, absorption and scattering of light to intermolecular interaction phenomena.

Simulation of Surface Delamination, effect on behaviour due to changes in relaxation radius and thickness

Simon Scarle, C. Ewels and M. I. Heggie
CPES, University of Sussex, Falmer, Brighton, BN1 9QJ.

Starting with a ball-and-springs approach, we present simulations of thin film or surface delamination behaviour. This work looks at the phase diagram of different delaminate structures produced by allowing only local relaxation out to fixed radii about failing bonds, and changes in film thickness.



Sequence specific binding of sodium ions to B-DNA studied by molecular dynamics simulation

Francesca Mocci and Giuseppe Saba
Dipartimento di scienze chimiche, Cittadella Universitaria di Monserrato, SS 554, bivio per Sestu, Cagliari Italy

Molecular dynamics simulations have been employed to probe the sequence-specific binding of sodium ions to the minor groove of B-DNA of three A*T rich oligomers having identical composition but different order of the base pairs : C(AT)4G, CA4T4G and CT4A4G. Recent experimental investigations, either in crystals or in solution, have shown that monovalent cations bind to DNA in a sequence specific mode, preferentially in the narrow minor groove regions of uninterrupted sequences of four or more adenines (A-tracts), replacing a water molecule of the ordered hydration structure, the hydration spine. Following these evidences it has been hypothesized that in A-tracts these events may be responsible for structural peculiarities such as a narrow minor groove and a curvature of the helix axis. The present simulations confirm a sequence specificity of the binding of sodium ions: Na+ intrusions in the first layer of hydration of the minor groove, with long residence times, up to ca. 3ns, are observed only in the minor groove of A-tracts but not in the alternating sequence. The effects of these intrusions on the structure of DNA depend on the ion coordination; when the ion replaces a water molecule of the spine, the minor groove becomes narrower. Ion intrusions may also disrupt the hydration spine modifying the oligomer structure to a large extent. However in no case intrusions were observed to locally bend the axis toward the minor groove. The simulations also show that ions may reside for long time periods in the second layer of hydration, particularly in the wider regions of the groove, often leading to an opening of the groove.

Simulation of the nematic-paranematic interface under shear

Guido Germano and Friederike Schmid
Theoretical Physics, University of Bielefeld, 33501 Bielefeld, Germany

Moving boundary conditions and a velocity profile unbiased thermostat [1] have been added to the domain decomposition molecular dynamics program GBmega, so that non equilibrium simulations of elongated molecules under shear flow can be performed on a massively parallel computer. This allows to simulate large scale systems including an interface between two coexisting phases like nematic and isotropic, that has recently been investigated by equilibrium simulations [2]. We present preliminary results on a system of 10^5 soft ellipsoids of elongation 15 in the nematic-paranematic coexistence region, as well as on smaller, homogeneous systems of the same particles, comparing them with theoretical predictions [3,4].

[1] D. J. Evans, G. P. Morriss, "Statistical Mechanics of Nonequilibrium Fluids" (Academic Press, London 1990).

[2] N. Akino, F. Schmid, M. P. Allen, Phys. Rev. E 63, 041706 (2001).

[3] P. D. Olmsted, C.-Y. D. Lu, Phys. Rev. E 60, 4397 (1999).

[4] P. D. Olmsted, Europhys. Lett. 48, 339 (1999).



Computer simulations of a side chain liquid crystal polymer

Lorna M. Stimson and Mark R. Wilson
Department of Chemistry, University of Durham, South Road, Durham DH1 3LE, U. K.

This study involves molecular dynamics simulations in the NpT ensemble of a model 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 microphases 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.



Computer Simulation Studies of a Liquid Crystal Dendrimer

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 dendritic liquid crystal systems. These multi-scale studies of a dendritic liquid crystal system model the third generation of a carbosilane dendrimer.

Firstly the system has 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 of the moment of inertia spheroid. We find that as a function of mean field strength the molecules change shape from spherical to rod-like.

Secondly a semi-atomistic 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.

Finally a coarse-grained representation of the LC dendrimer has been developed, incorporating structural information obtained from the earlier studies. This model has been studied in the bulk.



Calculation of Rotational Viscosity and Flexoelectric co-efficients of a nematic liquid crystal by atomistic computer simulation.

D.L. Cheung 1,2 S.J. Clark 2 Mark R. Wilson 1

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

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

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-efficient, γ1, and the flexoelectric co-efficients es and eb.

γ1 has been calculated using several methods: using the angular velocity correlation function of the nematic director Object1 and statistical mechanics methods based on the rotational diffusion coefficient 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. 99 9021 (1990)

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

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

Parameterisation and validation of a force field for liquid crystal molecules

D.L. Cheung 1,2 S.J. Clark 2 Mark R. Wilson 1

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

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

First principles density functional theory calculations have been carried out to determine the structure and conformational energies of a series of liquid crystal fragment molecules. These calculations employed a plane wave basis set for the valence electrons and the exchange-correlation interaction was evaluated at the Generalised Gradient Approximation level 1.

The calculations have been used to derive a molecular mechanics force field that describes a subset of liquid crystal molecules especially those used in display applications. In the parameterisation particular care was used to select functional groups that are not usually considered in other common force fields, such as laterally and terminally substituted biphenyls and partially conjugated units that are important in liquid crystal molecules.

The force field has been used to carry out molecular dynamics simulations 2 of the bulk phases of these liquid crystal fragments. Computed densities and heats of vapourisation are in good agreement with experimental data.

1. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais. Phys. Rev. B. 46 6672 (1992)

2. M. P. Allen and D. J. Tildesley. Computer Simulations of Liquids Oxford 1989



The prediction of helical twisting powers for liquid crystal chiral dopants

David J. Earl and Mark R. Wilson
Department of Chemistry, University of Durham, South Road, Durham DH1 3LE, U. K.

When a nematic liquid crystal is doped with a low concentration of chiral solute molecules a chiral nematic phase is induced in the system. The helical twisting power, $ \beta_{M}^{}$, of a chiral dopant is a measure of the twist the dopant is able to induce in a nematic phase. The ability to predict accurately the helical twisting power (HTP) of chiral dopant molecules is of fundamental theoretical interest in molecular chirality, as well as having useful technological applications in the design of new materials with high HTPs.

The work presented here uses Monte Carlo free energy perturbation calculations to determine HTPs. The method involves the growth of an atomistic model of a chiral dopant into a twisted nematic phase represented by soft repulsive spheroclinder1 (SRS) or Gay-Berne (GB) molecules in a sequence of simulations. The free energy change for this process can be determined over a full sequence of simulations. This same method is then employed for the enantiomer of the original molecule and $ \beta_{M}^{}$ can be calculated from the free energy difference for the growth of the chiral dopant molecule and its mirror image. Our results2, 3 show good agreement with experimental findings with the correct sign of twist being predicted in all cases. An alternative approach has also been studied where the twisted nematic liquid crystal is considered as a chiral potential rather than using explicit solvent molecules. Results for this method will also be presented here.

1. D. J. Earl, J. Ilnytskyi and M. R. Wilson, Molec. Phys. (2001), 99, 1719

2. M. J. Cook and M. R. Wilson, J. Chem. Phys, (2000), 112, 1560

3. D. J. Earl and M. R. Wilson, J. Mater. Chem., (2001), 11, 2672



Computer Simulation Studies of an Amphiphilic Polymer at an Air-Water Interface

Philip M. Anderson and Mark R. Wilson
Department of Chemistry, University of Durham, South Road, Durham DH1 3LE, U. K.

1,2-Dimethoxyethane (DME) contains the two major dihedral interactions present in the polymer poly(ethylene oxide) (PEO) (C-O-C-C and O-C-C-O). Consequently, its conformational energies have been extensively studied experimentally and theoretically in both the gas and liquid phases. This poster describes a combined quantum mechanical ab initio and classical molecular dynamics study to calculate the conformations of DME in the gas and liquid and to extract a force field to model PEO. The latter has been used to carry out molecular dynamics studies of an amphiphilic polymer at a water-air interface.

In this work, the ten major conformations (ttt, tgt, ttg, tgg, tgg', ggg, ggg', gg'g, gtg and gtg') and two major barriers (tct and ttc) in DME, have been subject to ab initio geometry optimisation, and energy evaluation at the MP2 and B3LYP levels of theory. These optimised energies were then used to parameterise a force field for DME, which was refined by carrying out molecular dynamics calculations of DME in the liquid phase, in an attempt to obtain good agreement with experimental data.

The force field was subsequently used in an atomistic molecular dynamics simulation of an amphiphilic graft-copolymer consisting of a polynorbornene backbone with PEO sidechains, sitting at a TIP4P water/air interface. At low surface coverage, in agreement with experimental observations, the polymer spreads out at the interface, but a small degree of penetration of the PEO chains into the bulk water is observed.

Density profiles from these calculations are being used to generate neutron reflectivity curves for this system, to compare with experimental neutron reflectivity studies.

Soft spheroid potentials for complex fluid simulations

Mike Allen 1 Guido Germano 2

  1. (University of Warwick)

  2. (Universitaet Bielefeld)

We discuss the use, in molecular simulation, of a class of soft-spheroid interaction potentials, closely related to the Gay-Berne potential and the spherical Gaussian potential which are commonly used in simplified liquid crystal and polymer simulations. We present a compact method of deriving forces and torques from the potential, re-deriving expressions first proposed by W Smith and K Singer, but not published in the primary literature by them. By applying thermodynamic perturbation theory, we show how these potentials may be used to simulate chiral liquid crystal molecules, or deformed shapes which exhibit flexoelectric behaviour.



A combined use of cluster approximations and exact information in the lattice-gas model for phase diagrams calculations

Yu.K.Tovbin, A.B.Rabinovich, D.V.Yeremich
Karpov Institute of Physical Chemistry, 10, ul. Vorontsovo Pole, 105064 Moscow, E-mail: tovbin@cc.nifhi.ac.ru

The use of numerical methods (Monte Carlo and molecular dynamics) to calculate the kinetics and dynamic characteristics of many processes in condensed phases under nonisothermal conditions is difficult, because this involves repeat determination of the concentrations of coexisting dense and rarefied molecule phases and the types of phase states, i.e., a knowledge of the entire phase diagram for a system. Such calculations are extremely time-consuming; therefore, it is necessary to develop methods that would provide calculation of phase diagrams and all equilibrium characteristics within a time typical of approximate methods, but would have accuracy close to that characteristic of the Monte Carlo and molecular dynamics methods. Therefore, it is reasonable to combine these two groups of methods [1]. The essence of the combine using the cluster and exact calculations is that the cluster solutions are fitted to exact ones by a variation of the calibration function parameters with simultaneously fulfilling the equations which describe the position of the critical point: dlnP/dq = d2lnP/dq2 = 0.

The Bether's approximation (or quasichemical approximation (QCA)) was considered as the cluster solution. (An improvements of the QCA with a help of the calibration function were discussed early in papers [2,3(14-16)]. The QCA is better than mean field approximation as well as it allows calculating the phase diagrams very quickly. Such combined modification gives the increasing accuracy in area of critical point within a time typical of approximate methods.

The combined approach was tested on square lattice by comparisons with exact Onsager solution [4,5] for coexisting curve and with Monte Carlo results for c(2x2) ordering system [6]. The approach was employed for the description of phase diagrams for two-dimension lattice at adsorption systems on homogeneous and heterogeneous surfaces (the phase diagrams for argon adsorption on rutile and muscovite), and also for three- dimension system at adsorption in porous adsorbents (argon adsorption in narrow slit-like carbon pores) [7]. The exact information on areas of critical points in homogeneous lattices was used. The difference between homogeneous and heterogeneous lattices was described by the cluster approximation, which gives enough accuracy to extrapolate of the exact solution on homogeneous lattice for heterogeneous lattices. A modern state of simulations of the phase diagrams in porous systems is discussed.

The using of the calibration functions permits to improve the description of phase diagrams and various adsorption characteristics in vicinity of critical points without increasing of a computation time.

Acknowledgment. The work was carried out at financial support of INTAS (project 00-505).

References

1. Yu.K.Tovbin, Russ. J. Phys. Chem., 72 (1998) 2280.

2. I.P.Prigogine, L.Mathot-Sarolea and L.van Hove, Trans. Farad. Soc., 48 (1952) 485.

3. P.Cavalloti, Surface Sci., 83 (1979) 325.

4. L.Onsager, Phys Rev., 65 (1944) 117.

5. R.J.Baxter, Exactly Solved Models in Statistical Mechanics, London: Academic Press, 1982.

6. K. Binder and D.P.Landau, Surface Sci., 108 (1981) 503.

7. Yu.K.Tovbin, A.B.Rabinovich and E.V.Votyakov, Russ. Chem. Bull., 51 (2002) No.9, accepted for a publication



Mechanism of palladium-diphosphine catalysed copolymerisation of olefins and carbon monoxide

Christopher Neville-Smith
Chemistry, Durham University (Trevelyan College), Supervised by Mel Kilner

Palladium diphosphine catalysts are known to copolymerise olefins and carbon monoxide under suitable conditions to form copolymers of the form (-CH2CH2CO-)n. These copolymers are created through perfect alternation of addition and insertion of carbon monoxide and the olefin. There is considerable interest in the mechanism of the reaction, in particular the reasons for the lack of double insertion of olefins or carbonyls. A number of different reasons for alternating insertion have been suggested, but it is unclear how important each of these factors are.

Previous work by Zeigler and Margl, and Svensson et. al., have modelled the CO and olefin insertion steps in some detail, but did not consider the mechanism of the two addition steps in any depth. This project aimed to optimise the mechanism of the entire propagation cycle, with scope for later adding the initiation and termination steps.

So far, the key findings have been as follows:

Theory and simulation of the nematic anchoring coefficient

Denis Andrienko 1 and Michael P. Allen 2

  1. Max-Planck Institut fuer Polymerforschung, Ackermannweg 10, 55128 Mainz, Germany
  2. Department of Physics and Centre for Scientific Computing, University of Warwick, CV4 7AL, United Kingdom

Combining molecular simulation, Onsager theory and the elastic description of nematic liquid crystals, we study the dependence of the nematic liquid crystal elastic constants and the zenithal surface anchoring coefficient on the value of the bulk order parameter