Abstract

Stochastic classical molecular dynamics coupled to functional density theory: applications to large molecular systems

K. C. Mundim

Institute of Physics, Federal University of Bahia,

Salvador, Ba, Brazil

and

D. E. Ellis

Dept. of Chemistry and Materials Research Center

Northwestern University, Evanston IL 60208, USA

Received 07 December, 1998

I Introduction

In complex materials and biophysics problems the number of degrees of freedom in nuclear and electronic coordinates is currently too large for effective treatment by purely first principles computation. Alternative techniques which involve a hybrid mix of classical and quantum methodologies can provide a powerful tool for analysis of structure, bonding, mechanical, electrical, and spectroscopic properties. In this report we describe an implementation which has been evolved to deal specifically with protein folding, pharmacological molecule docking, impurities, defects, interfaces, grain boundaries in crystals and related problems of 'real' solids.

As in any evolving scheme, there is still much room for improvement; however the guiding principles are simple: to obtain a local, chemically intuitive description of complex systems, which can be extended in a systematic way to the nanometer size scale. The computational scheme is illustrated in Fig. 1, showing classical Molecular Dynamics (MD), Monte Carlo (MC-GSA) stochastic sampling, and Discrete Variational (DV) Density Functional (DF) quantum mechanics procedures coupled together. We next briefly describe each of the component procedures.

II Classical Dynamics, Stochastic Methods, and Quantum Clusters

II.1 Classical Methodology: Molecular Dynamics

The idea of molecular modeling is an attempt to describe the quantum chemical bonds in terms of a classical force field in Newton's equations. In a Molecular Mechanics (MM) model the atomic bonds are represented by springs joining the atoms; i.e., the molecule is assumed to be a collection of masses and springs. In this case is necessary find a set of parameters, for this force field, that fit the quantum atomic interactions with sufficient accuracy. Ideally one hopes that such refinements will eventually lead toward a unified computational model that can successfully mimic observed molecular properties. Academic research efforts and the pharmaceutical industry interest in developing new compounds in biological molecular systems have stimulated the appearance of different computational codes based on classical force fields.

The first molecular dynamics applications was performed by Alder and Wainwright, [1] who used a perfectly elastic hard sphere model to represent the atomic interactions. One of the most widely used force field ones is called MM1 (Molecular Mechanics), proposed by Allinger [2]. At the moment there are an uncountable number of MM computational codes using a classical force field. Each one of these uses a particular force field to describe different molecular properties and to fit some experimental results.

Reasons which justify the increasing use of MM can be listed, as for example:

- Possibly the chief reason is the relatively short computational time, which for MM methods increases as M², where M is the number of atoms in the molecule. In contrast, the use of the ab-initio quantum methods in such molecular systems is computationally impracticable because the computer time to evaluate the inter-electronic repulsion integrals increases as N⁴ (or more rapidly, when correlation is taken into account), where N is the number of basis orbitals. Normally, there are at least several basis functions per atomic orbital shell- 1s, 2s, 2p, etc..

- The MM method and its results are conceptually easier to understand than quantum mechanical (QM) methods.

- In MM, it is very simple to introduce time evolution;

- In MM, it is possible to introduce the temperature as an external perturbation, and trace its effects.

Two important questions arise which are unfavorable to the MM methods: First, there do not exist well-defined rules to evaluate the force constants; and second, in order to choose the best force field is necessary to have a priori some considerable knowledge about the molecular system. In molecular mechanics the atom is represented by a spherical body with a particle mass equal, in general, to the respective atomic mass. In today's molecular mechanics, several force field models have been proposed. In general the molecular energy (potential function) related with the classical force field can be expressed by the following sum [3],[4];

Where, in summary, the first term of eq. (1) V_H represent the energy necessary to stretch or compress the atomic bond; Vb is the potential contribution that represents the bond-angle bending interaction; V_tor-i, V_tor-p are the torsional contributions that represent the harmonic dihedral bending interaction and the sinusoidal dihedral torsion interactions. The last terms V_C and V_vdW represent the non-bonded interactions such as Coulombic repulsion, hydrogen bonding and van der Waals interactions. Each of these potential energy functions represents a molecular deformation from a reference geometry. The interactions given in equation (1) are available in our Molecular Mechanics code, and are typical of other currently available codes. To study molecular properties in solid systems, additional potentials have been introduced, including Born-type exponential repulsion and Morse potentials.

Molecular dynamics (MD) calculations consist of analyzing the evolution of the molecular system with time. In this case the atoms are continuously moving, the bonds are vibrating, the angles are bending and the whole molecule is rotating. In MD, successive configurations of the system are generated by integrating Newton's equations of motion. The result is a trajectory that specifies how the positions and velocities of the atoms vary with time. The trajectory is obtained by solving the following differential equations;

where is the force over atom i along the coordinate and m_i is the atomic mass.

Equation (2) is solved by a difference finite technique with continuous potential models (1), which are generally assumed to be pairwise additive. The essential idea is that the integration is broken down into many small stages, each separated in time by a fixed increment dt. For each step the total force on each atom is calculated as the vector sum of the interactions with some or all atoms belonging to the molecular system. The force is calculated by taking the gradient of the potential function represented by eq. (1). Using the forces and accelerations, all atomic positions and velocities are determined using Newton's equation at the time t + dt. In the same procedure the forces and the new velocities are used to calculate the new atomic position at time t + dt, and so on.

Once the force on each atom has been calculated for time t, the position of each atom at some later time t + dt is then given by

where . Here t = ndt, and dt ~ 0.001 pico seconds. The initial coordinates may be taken from experimental data or may be randomized; the initial velocities are derived from the Maxwell-Boltzmann distribution, in general, and temperature T is controlled by a dynamic scaling procedure. A cooling or annealing procedure of modifying T provides one way to approach stationary points, and possibly the equilibrium state of the system. Once the potentials are calculated for each new geometry, the new forces and positions of each atom can again found. This procedure is repeated until convergence in the energy is reached.

There are different algorithms or numerical methods available to solve differential equation (2). All methods are based on finite differences and solve the equation step by step in time. Often the step size is taken constant, but adaptive methods have sometimes been found useful.

The simplest and most straightforward (but unfortunately not sufficiently accurate) approach is to use the Taylor expansion based on equation (2).

Another method sometimes used (but computationally very expensive) is the Taylor predictor;

This algorithm contains no explicit velocities, and we see that also the third derivative terms can be made to cancel. In this case the velocities can be approximated by

Another very well known algorithm is the Verlet leap-frog scheme given as;

One consideration for molecular dynamics that immediately invalidates a number of these methods is the cost of calculation of the force, or gradient of the potential. Computation of the force is extremely laborious compared to any manipulation involved in updating the variables to take one step forward in time. This means that any method that involves more than one force evaluation per step cannot be a good choice. This rules out the well-known Runge-Kutta method and its variants, that require up to four force evaluations per step, so the simple Verlet scheme is generally preferred.

The MD procedure can be briefly summarized in the following steps:

i- Guess the initial molecular geometry, the initial atomic velocities and environment temperature.

ii- Calculate the force, as the gradient of the potential (1).

iii- Calculate the new atomic position or molecular geometry using one of the algorithms to integrate Newton's equations.

iv- Redefine the atomic velocities using the virial theorem

Where k_B is the Boltzmann's constant, T_o is the initial or "bath'' temperature and T is the local temperature. v_old and v_new are the velocities of two consecutive steps in the MD procedure.

v- Save intermediate information for statistics and go back to step (i) until the stabilization of the molecular energy is reached. At the end of the MD procedure the total energy (kinetic and potential contribution) must be constant.

This methodology, though pedagogically correct, allows the system to get "trapped'' in local energy minima. There is no guarantee that this geometry corresponds to the global minimum energy. So the usual MD procedure may not be the best way to search for the global minimum point on the complex energy hypersurface.

In this context, if the molecular hypersurface energy is non-convex, stochastic molecular dynamics offers a more efficient way of finding both local minima as well as the global one. Thus we next discuss the more important concepts needed to implement a stochastic molecular dynamics procedure.

II.2 Stochastic Dynamics: Generalized Simulated Annealing

Simulated Annealing [5][6] methods have been applied successfully in the description of a variety of global extremization problems. Simulated Annealing methods have attracted significant attention due to their suitability for large scale optimization problems, especially for those in which a desired global minimum is hidden among many local minima. The basic aspect of the Simulated Annealing method is that it is analogous to thermodynamics, especially concerning the way that liquids freeze and crystallize, or that metals cool and anneal. At high temperatures, the molecules move freely with respect to one another. If the system is cooled slowly, thermal mobility is lost. The atoms are often able to line themselves up and assume a molecular geometry that is in general a local equilibrium state. The simulated annealing procedure is actually more complicated than the combinatory one, since the familiar problem of long, narrow potential valleys again asserts itself. Simulated annealing, as we will see, tries random steps, but in a long, narrow valley, almost all random steps are uphill. The amazing fact is that, for a slowly cooled system, nature is able to find this minimum energy state. So the essence of the process is slow cooling, allowing ample time for redistribution of the atoms as they lose their mobility. This is the technical definition of annealing, and it is essential for ensuring that the lowest energy state will be achieved.

The first nontrivial solution along this line was provided by Kirkpatrick et al. [5],[6] for classical systems, and also extended by Ceperley et al. [7] to quantum systems. It strictly follows the quasi-equilibrium Boltzmann-Gibbs statistics using a Gaussian visiting distribution, and is sometimes referred to as Classical Simulated Annealing (CSA) or the Boltzmann machine. The next interesting step in this subject was Szu's proposal [8] to use a Cauchy-Lorentz visiting distribution, instead of a Gaussian distribution. This algorithm is referred to as the Fast Simulated Annealing (FSA) or Cauchy machine.

On the other hand, a Generalized Simulated Annealing (GSA) approach which closely follows the recent Tsallis statistics [9],[10] has been proposed [11],[12],[13],[14]; GSA includes both the FSA and CSA procedures as special cases. We have implemented the GSA algorithm as a method to calculate the minimum energies of conformational geometries for different molecular structures. This technique can be applied in either quantum [12] or classical [13] methods. The GSA method is based on the correlation between the minimization of a cost function (conformational energy) and the geometries randomly obtained through a slow cooling. In this technique, an artificial temperature is introduced and the system is gradually cooled in complete analogy with the well known annealing technique, frequently used in metallurgy when a molten metal reaches its crystalline state (global minimum of the thermodynamic energy). In our case the temperature is intended as an external noise, which acts as a convenient stochastic source for eventual detrapping from local minima. Near the end of the process the system hopefully is within the attractive basin of the global minimum. The challenge is to cool the system as fast as possible and still have the guarantee that no irreversible trapping at any local minimum has occurred. More precisely, we search for the quickest annealing (approaching a quenching) which maintains the probability of finishing within the global minimum equal to one.

The procedure to search the minima (global and local) or to map the energy hypersurface consists in comparing the conformational energy for two consecutive random geometries x_t+1 and x_t obtained from the GSA routine. x_t is a N- dimensional vector that contains all atomic coordinates (N) to be optimized. The geometries, for two consecutive steps, are related by

where Dx_t is a random perturbation on the atomic position.

To generate the random vector Dx_t the present GSA routine uses an extension of the procedure given in Ref. [13]. We have calculated Dx = g^-1(w) using a numerical integration of the visiting distribution probability g_qv(x). Where w is a random vector [0,1] obtained from an equiprobability distribution and g^-1 is the inverse of the integral of g_qv(x) given by

Mathematical details of the structure of the distribution function g and its inverse g-1 are given in the reference [13]. In summary, the complete algorithm for mapping and searching for the global minimum of the energy is:

(i) Fix the parameters (q_A; q_v). We note that (q_A;q_v) = (1;1) and (1;2) respectively correspond to the Boltzmann and Cauchy machines. Start at t = 1, with arbitrary internal coordinates and high enough value for visiting temperature T_qv(1) and cool as follows:

where t is the discrete time corresponding to steps of computer iteration.

(ii) Next, randomly generate the new atomic coordinate x_t+1 from x_t as given by the visiting distribution probability gqv as follows:

For sufficiently long time simulations this procedure assures that the system can both escape from any local minimum and can explore the entire energy hypersurface. This equation is used in the GSA routine and differs from the general proposal given in [11]; instead we build a minimization vector using (8).

(iii) Then calculate the conformational energy E(x_t+1) from the new molecular geometry using the classical force field [3]. The new energy value will be accepted according to the following rules:

if E(x_t+1) £ E(x_t), replace x_t+1 by x_t; if E(x_t+1) > E(x_t), run a random number r Î [0,1];

if r > P_qA (acceptance probability) retain x_t; otherwise, replace x_t by x_t+1.

The acceptance probability is given by

(iv)- Calculate the new temperature T_qv(t) using Eq. (12) and go back to (ii) until the convergenceof E(x_t) is reached within the desired precision.

In order to make clear the procedure to construct the presently used computational code, we present the following flowchart as Fig. 2:

II.3 Quantum Methodology: Density Functional Embedded Cluster Scheme

The Density Functional (DF) theory is a standard tool for describing electronic structures, which is first-principles in concept and flexible in execution. DF theory has become a major tool for analyzing electronic structure of molecules and solids, with high intrinsic accuracy and reasonable computational effort. For extended systems with no periodicity to exploit, a cluster approach offers many advantages, providing that a reasonable embedding scheme is used. Definition of 'reasonable embedding' depends upon the nature of the system- while a point charge environment is acceptable for a highly ionic compound, it is quite poor in treating localized covalent interactions where there is hybridization of orbitals with neighboring atoms. The embedded cluster density functional (ECDF) scheme permits analysis of wavefunctions and their derived properties in a restricted volume, with interactions to the extended host included by effective potentials and boundary conditions. The present implementation of the ECDF scheme employs 'guard atoms' on the surface of the cluster to saturate interior bond structures, and synthesis of total charge and spin densities including the environment, to determine effective cluster potentials. Scanning algorithms by which the 'view window' of the cluster is moved over the volume of interest, permit treatment of extended systems. A key concept in obtaining self-consistency of the multiple overlapping views of a complex system like a heterogeneous interface, is that of equilibration of the chemical potential, or Fermi energy EF, among the component clusters. This methodology has been developed in several forms over the years, and successfully applied to a variety of molecules and extended solids, including metals, semiconductors, and insulators; see for example [15].

In outline, the following main steps are performed in ECDF calculations:

i) Charge and spin densities of the extended system are constructed by summing component atom-like densities for the cluster and host:

ii) LCAO-type basis functions are generated by numerical solution of DF equations for atoms/ions in a potential well, determined from the previous self- consistent-field iterations; the basis is thus iterated as part of the overall optimization procedure.

iii) Form matrix elements of the effective Hamiltonian

and overlap operators over the AO basis c_j , possibly transformed into symmetry orbitals, and orthogonalized against a frozen core. Here t, V_C , and V_xC are kinetic energy, nuclear and electronic Coulomb potential, and exchange-correlation potential respectively. The index s refers to spin orientation. The DF potential is formed from the densities, and boundary conditions are applied to force localization of solutions to the cluster.

iv) Solve the Schrödinger (nonrelativistic) or Dirac (relativistic) equation, to obtain a variational expansion of the one-electron wavefunctions as

v) Project the cluster densities onto a localized 'effective atom' expansion as

Here f_n are Fermi-Dirac occupation numbers and u enumerates atoms within the cluster.

vi) Analyze the interior 'seed atom' region of each cluster to extract component atomic densities. New atomic densities for the next iteration are obtained after equilibration of the system and application of constraints such as electroneutrality of crystalline unit cells. Equilibrate the various clusters selected to span the region of interest, either by matching Fermi energies or spectral features, and iterate steps i-vi until adequate convergence of charge and spin distributions, spectral features, and total energy is obtained.

vii) Extract densities of states, optical and X-ray absorption coefficients, electron energy loss spectra (EELS), cohesive energies, and other desired properties. Compare charge distribution and cohesive energies with data assumed in construction of MD/GSA potentials. Update MD/GSA potentials as required.

III Examples of Applications

Each of the components of the hybrid scheme: MD, MC/GSA, and ECDF have been developed independently and has a considerable history of successful applications. The linkage between these components, with the capability of feeding information between the components provides us with a powerful new tool for studying materials properties. An area of intense development in our laboratories is design and construction of a graphical user interface which will permit the use of this system of programs by non specialists, and particularly by experimentalists wishing to analyze new materials. In the remainder of this report, we give four illustrative examples related to current problems in materials design and optimization.

Beyond the general procedures outlined above, successful applications depend upon some degree of experience and 'art'. For example; when is it most efficient to use MD and when is it better to use MC/GSA to determine the minimum-energy configuration of a system as complex as a ceramic grain boundary containing impurities? In practical terms, we observe an initial rapid convergence with the MC/GSA approach, followed by a long 'tail' of sampling in which only small incremental improvements are found. Thus an optimal scheme involves switching from GSA to MD-gradient based procedures at a rather well defined point in the simulation, when an accurate energy minimum is desired. Further details will be given elsewhere; here we only wish to comment briefly on another important practical aspect: the use of boundary conditions and tight constraints to permit efficient sampling of regions of interest. In traditional MD/MC one often tries to work with the largest possible simulation volume, to minimize boundary effects, and distortions due to imposed periodicity. Here, we deliberately choose quite restricted simulation volumes, with tight boundary conditions, in order to focus upon system response to perturbations related to defects, surfaces, and interfaces. Of course, too tight control completely prejudices the outcome of simulation, so convergence and stability tests have to be made as in any local scheme.

During the dynamic process the MD/MC procedures may call the SCF-DV procedure to redefine the atomic charges. In case of ionic interactions, e.g. molecules docking with ionic lattices, this offers the opportunity to determine polarization and dynamic charges 'on the fly' as a function of geometry. In general, the DV quantum results produced can be used to adjust the classical force field (See flowchart in Fig.1); however, the interaction among various potential parameters has to be considered. For example, in a typical Born-type ionic interaction with point-charge long range terms, exponential repulsion, and inverse power short range attraction, modification of point charge values also implies some modification of short-range parameters to maintain correct bond lengths and cohesive energies.

III.1 Applications Using a Coupled DVM-GSA Method

III.1.1 A Stacked Linear Porphyrin Chain

The chemical bonding of metal porphyrins, porphyrazines, and phthalocyanines is of importance in understanding biophysical and catalytic processes. The crystalline materials like copper phthalocyanine (CuPc) form stacked chains, and when partially oxidized by iodine (CuPcI) for example, become interesting 'molecular wire' conductors. Doped systems like Cu(1-x)Ni(x)PcI show quasi-one dimensional magnetic interactions of considerable theoretical interest. Component molecules (monomers) like CuPc are quasiplanar and small enough to be treated by standard quantum chemical approaches; however, understanding the electronic coupling between monomers requires a DF approach, and understanding the dynamics of coupling requires a MD/MC classical treatment. An ECDF/MD/MC study of a CuPc stack has been recently carried out [16]; here we wish to give a few highlights of the results.

Fig.3 represents the electrostatic potential in a 3D-mapped isosurface of a porphyrin chain. The molecular equilibrium conformation was obtained using a GSA technique, coupled to a classical force field. We have carried out semi- empirical Hamiltonian (PM3) calculations to draw the electrostatic potential isosurface. One of the significant conclusions drawn from this work is that standard force fields are capable of reproducing both monomer and intermolecular interactions with considerable accuracy. Thus, while the underlying electronic interactions (porphyrinic-ring pi versus metal-d; sigma bonding in-plane versus pi-interaction between molecules) are complex and require detailed analysis, a simple classical parametrization captures the main structural results. This observation helps to justify and motivate the mixed classical/quantum methodology.

Figure 3. Electrostatic potential in a 3D isosurface.

Figure 3a. Linear stacking of a Porphyrin Chain.

III.1.2 Carbon Interstitials in Copper

In metallic conductors like copper, small particle precipitates are found to improve hardness and wear resistance of contacts. For high temperature applications, matrix-embedded carbon fibers have also been considered, to provide stiffness and prevent ductile failure, and may be useful to reduce sliding friction [17]. An understanding of C diffusion in Cu, and its behavior at the Cu/C interface is fundamental to developing successful strategies for improving fiber-matrix composite performance. Wettability-enhancing efforts based upon alloying elements benefit by connection between the fundamentals of the wetting of solids by liquid metals and the practice of the preparation of metal-matrix composites (MMC) [18]. In [19] it was shown that the interfacial energy in a metal-carbon fiber system is modified by the interfacial adsorption of the alloying elements which are added to the metal matrix.

High-resolution scanning electron microscopy (HRSEM) shows the existence of a 50nm solid solution zone at the copper-carbon interface, after annealing. After heat treatment of one hour at 1000 ^oC the interface faded and we can see a copper-carbon solid solution region of width 50 nm. Formation of a very dilute Cu-C interstitial solid solution at the interface may be diffusion- controlled. To predict the diffusion behavior of carbon atoms we have to study the structure and interatomic interactions in such alloys. Phenomenological studies of these phenomena have been augmented by quasi-continuum models [20] which reflect some aspects of the underlying electronic structure. Now the analysis is extended to the atomic scale, using a combination of GSA/MD atomistic simulations and embedded cluster DF schemes.

Our MD simulations [20] are here deliberately constrained to small displacements of the Cu host atoms, for which the framework of elastic Hooke's- law Cu-Cu interactions is quite adequate. As follows from simple pseudopotential calculations, the Cu-C interaction is repulsive at least to several coordination shells surrounding the carbon atom. The repulsive part of interaction may be fitted by any functional form, and we normally choose Lennard-Jones (LJ) parameters to reproduce the experimentally measured properties. Unfortunately, experimental data on dilute Cu-C solid solutions are absent, so we use the results of the previously mentioned pseudopotential studies.

The expression used for the classical system energy is thus

The last term, E₀, represents the volume-dependent part of the energy and includes information about electronic density redistribution, which is gathered from the embedded cluster density functional scheme, or may be obtained from other electronic structure calculations. The parameter values are: R₀=2.475Å , K=6.920 eV/Å ², E₀=-0.875 eV, C6=41.548 eV/Å ⁶, C12= 2989.105 eV/Å ¹².

Among the major conclusions of this work we find the following:

i. Carbon strongly prefers the octahedral interstitial site over tetrahedral sites, as predicted from earlier works. We have determined the relative site energies and the (rather long range) relaxation of the Cu host around the impurity.

ii. Surface sites (incomplete octahedral interstitials) are at lower energy than volume sites, consistent with the extremely low solubility of C in Cu.

iii. Carbon dimers and hints of clustering appear at higher impurity densities.

iv. The graphite-Cu interface region appears to be disordered to a depth of several Cu atomic layers.

III.1.3 Scheelite, a Candidate Host for Matrix-Fiber Composites

Calcium tungstate, CaWO₄, also known as scheelite, is a prototype for a large class of naturally occurring minerals of the same structure [22]. The WO₄ tungstate group exhibits strong bonding of mixed ionic and covalent character of great rigidity; as a result the scheelite lattice has a rather anisotropic response to pressure, with compression along the c-axis (see Fig. 5) being ~ 1.2 times that in the a-b plane. The W-O bond length is 1.75Å ; in contrast, the Ca-O bond length is 2.41Å and the Ca environment may be described best as an eight-fold coordinated cage, formed from oxygens of eight different surrounding WO₄ tetrahedra.

The contrasting bond structures are the underlying reason for the apparently contradictory high melting temperature and low fracture strength of the material. These are precisely the characteristics which make scheelite interesting as a host matrix for ceramic matrix-fiber composites (CFC). A useful CFC should be air stable to T > 1300 ^oC, and under stress, cracking should take place in a controlled manner along the matrix-fiber interface, allowing the fibers to pull out and minimize structural damage. A typical fiber candidate would be a-Al₂O₃ and and so the surface structures of scheelite and alumina, and their interfaces are of primary importance. Various crystallographic surfaces of alumina have been characterized by experiment and theory; however, very little is known about the low energy surface structures of CaWO₄. Simple considerations suggest that oxygen-terminated surfaces which preserve the WO₄ units should be of relatively low energy; thus we have carried out GSA modelling of the (001) oxygen terminated scheelite surface. A portion of the hemispherical (001) surface volume model is shown in Fig. 6; rigid embedding atoms included in the force field are not shown.

The potential parameters used are given in Table 1. Note that the strong W-O 'bond spring' and angular spring constants are chosen to reproduce the observed rigidity of the WO₄ structural unit. The expected result is that the greatest response to crystal cleavage will be some rearrangement of the Ca ions with respect to the tungstate network. The surface relaxation energy is calculated to be 17 meV/Å ² with respect to the cleaved unrelaxed crystal, with tungstate angles held fixed, and 20 meV/Å ² when flexure is permitted.

The (001) oxygen-terminated cleavage plane, with its intact WO₄ groups, appears as the most obvious low energy surface. There are several other cleavage planes parallel to the c-axis which may be of relatively low energy. Interestingly, simulations of these surfaces by MD techniques suggest that cleavage along such places will result in spontaneous rupture along the (001) plane [23].

We have carried out ECDF calculations on the bulk and surface environments to determine major features of the electronic structure. The simplest analysis of the results can be made in terms of Mulliken atomic orbital populations and densities of states (DOS) diagrams. In Table 2 we give the self-consistent orbital populations and charges for bulk and surface regions, noting that the formal valencies are C^a+2, w⁺⁶, O^-2. Consistent with our qualitative discussion above, we see that calcium is near its nominal ionic value, indicating the relatively weak and long-range nature of its interactions. In contrast, the effective charge of ~ +3 on tungsten emphasizes the covalent charge sharing in the (WO₄)²- group. As is generally found in metal oxides, the oxygen 2p band, while formally fully occupied, has in fact 0.5e vacancies per atom, due to charge transfer from the metals. The bonding-antibonding gap in energy levels of the crystal-embedded (WO₄)²- complex is found to be ~ 5.0eV, indicative of its great stability. The occupied valence band region shows W 5d, 6sp structures spanning ~ 6eV and forming two distinct sub-bands with the O 2p, which can be qualitatively labeled as p- and s-type or equivalently denoted as e- and t- symmetry of the T_d local point group.

In order to begin to explore the oxide-oxide interfaces critical for understanding matrix-fiber interactions, we have generated a preliminary model for CaWO₄ (001): MgO(001). To induce steps and irregularities, the two surfaces have been made to intersect at an angle of ~ 25 degress. The volume between the two perfect bulk crystals was 'regrown' using the GSA scheme with variable atomic composition and positions. ECDF studies were made of selected interface regions to examine rearrangement of electron densities associated with changes in coordination and bond lengths. A preliminary report on this work has been presented [24]; a detailed analysis of cohesion and bonding effects will be given elsewhere. A general conclusion that can be drawn from the surface and interface studies is that reduction in cation coordination and concomitant reduction in average metal-oxygen bond lengths leads to reduced ionic charges at solid-vacuum and oxide-oxide interfaces. In the case of scheelite and related tetrahedrally bonded materials like monazite, LaPO₄, the robust tungstate (phosphate) groups resist deformation and reduction processes, leading to special fracture behavior and high temperature chemical stability.

Finally, we mention results of recent simulations on the scheelite (100): a- Al₂O₃ (0001) interface [23]. Rather large unit cells were used in slab geometry with two-dimensional periodicity, using a gradient-search MD scheme. It appears that relaxation is essentially localized to a few atomic layers on either side of the interface, due to the considerable rigidity of the Al-O and W-O bonds. The main relaxation mechanism is the rotation of the WO₄ groups to accommodate interface strain. ECDF calculations have been undertaken to verify the local electronic structure and to reconcile this structure with the force field model.

III.2 Molecular Modeling in New Drug Development

The development of a new pharmaceutical is a long and expensive process. A new compound must not only produce the desired response with minimal side- effects, but must be demonstrably better than existing therapies, and produced at a reasonable cost. Finding novel compounds with desirable properties can be a difficult problem, and one where molecular modeling has much to offer. The development of molecular modeling techniques is contributing to the understanding of biological processes at the atomic-molecular level and to the proposal of new molecular structures with high biological efficiency. The action of many drugs, hormones and membrane proteins is dependent on the spatial conformation adopted by the molecules in the inhomogeneous environment formed by the cell membrane and its neighborhood.

The modeling of biological molecules has led to great progress in the design of new molecular structures with desired specific properties. Such products could be obtained synthetically or by genetic engineering procedures. Molecular modeling can lead to the understanding of the dynamic circumstances compatible with experimental observations, permitting a formal definition of the conditions that would be expected to produce the observed behavior and therefore to infer patterns of behavior for situations of interest.

III.2.1 Conventional Molecular Dynamics Procedure Applied to Biomolecular Systems

Our computational code, based on the classical force field presented in section (2A), has an additional capability of modeling molecules of biological interest in aqueous and apolar media, as well as at their interface. The software establishes two media with different continuous dielectric constants, separated by a cytoplasm/membrane environment. The approach used to simulate the membrane interface is expressed by a discontinuity in the dielectric constant, taking into account the different electrical polarizability of the aqueous and the hydrocarbon phases. Then, the electrostatic interaction between non bonded atoms at this interface will be renormalized by the influence of the dielectric discontinuity. The polarization field produced at the surface of discontinuity by a point charge can be calculated by the method of images. A fictitious charge is placed in the opposite phase: the distance and the value of the image charge are fixed by taking the appropriate electrical boundary conditions at the surface.

When two charges i and j are on the same side of the interface, the potential on i due to j will have a contribution due to the image of j (i.e., j'), so that the Coulomb potential on i will be expressed by

where

and where r_ij and r^¢_ij are the distance between i and j, and i and j' respectively. X_s is the position of the surface between the two media. When two charges are at different sides of the interface the method of image gives the following result:

We have applied this MD approach to study to study a new mutant sequence of E. Coli. , modeling a 21 amino acid peptide, a mutant sequence from E. Coli maltoporian [3],[4]. An important aspect of peptide simulation is the great number of degrees of freedom for linear conformations, which can lead to several conformations with equivalent values of minimum energy. In this investigation we have studied a mutant signal sequence (D78r₁) of the LamB gene product, also called l receptor or maltoporin, a protein of the external membrane of E. coli. The mutant sequence contains 21 residues and shows a deletion of four residues in the hydrophobic region relative to the wild sequence. Such a deletion should abolish its capacity for helix formation and therefore its functionality. These effects were confirmed experimentally. However, 50 was re-established by replacing a Gly residue with a Cys residue at residue 13 of the mutant peptide. The relevant sequences are

Polarization effects on the peptide conformations have been investigated through the electrostatic charge image method as described above. A similar technique to simulate polarization effects in peptide conformations has been proposed by some of the authors in Ref. [25]. These authors carried out a classical force field simulation containing electrostatic image charge calculations to investigate the conformations of peptides characterized by different hydrophobicities at a water-membrane interface model. The interface is represented by a surface of discontinuity between two media with different dielectric constants, taking into account the difference between the polarizabilities of the aqueous medium and the hydrocarbon one.

III.2.2 Stochastic Molecular Dynamics Procedure Applied to Biomolecular Systems

The biological function of a protein or peptide is often intimately dependent upon the conformation(s) that the molecule can be adopt. X-ray crystallography and NMR are two methods used to provide detailed information about protein structures. Unfortunately, the rate at which new protein sequences are determined experimentally is much greater than that for determination of structures. There is thus considerable interest in theoretical methods for predicting the three-dimensional structure of a protein from its amino acid sequence; this is called the protein folding problem. Prediction of protein folding is a case where the usual MD is not a good theoretical method, since the real physical time of a folding process may be on the scale of minutes. To solve the time evolution of Newton's equation in MD it is necessary to discretize the time in steps of dt ~ 1 femtosecond, due to the intrinsic vibrational frequencies. This means that about 6x10¹⁶ MD steps would be needed to simulate the protein folding process, which is beyond computational possibility.

An alternative solution to this problem is to use stochastic simulation procedures such as GSA. Recently, Moret [26] has discussed the efficiency of stochastic simulation to study protein folding by the GSA method. The author simulated a-helix configuration of the ico-alanine polypeptide (Fig. 7) and found the necessary conditions for GSA parameters to obtain adequate computational efficiency in protein folding simulation. The GSA procedure with visiting index q_V =1,05 was shown to be more efficient than the usual simulated annealing method of Kirkpatrick (q_V = 2) and the fast simulated annealing method of Szu (q_V = 1) for calculation of an a-helix conformation equilibrium structure.

IV Conclusions:

A nearly optimum procedure to search for minimum-energy structures is to initially scan using the MC/GSA scheme. As the slow-convergence region is attained, we switch to the MD/annealing scheme. The algorithms are coupled together in such a way as to permit this strategy, in a unified and convenient manner.

To extract electronic structure information in 'interesting' nuclear configurations, and not only in the state of lowest potential energy, we have the capacity to activate the Embedded Cluster Density Functional components of the codes from within the running dynamics procedures. This permits generation of 'snapshots' of electronic states during a process of molecular interaction or docking of a molecule with a surface, wall of a molecular sieve, or macromolecule. An important application which helps in dynamically re-parametrizing the interatomic potentials, is the determination of atomic charges 'on the fly'.

With the concept of self-consistent embedding, we have gained the ability to resolve large-size structures, containing thousands of atoms, and a size scale of at least 10 nm, as we have illustrated in the examples given here. Studies in progress on oxide surfaces, molecule-zeolite interactions, and grain boundaries in metals will help to determine future evolution of the hybrid classical/quantum procedure. Although GSA/MD/ECDF is already a powerful tool for the study of complex systems, it is not optimized for processes which are essentially dynamic, such as diffusion of atoms and molecules along interfaces and surfaces. It will be important to develop the capability for treating the longer time scales implied by diffusion processes; it is encouraging that the protein folding problem (which has a long time scale) has already been shown to be amenable to the GSA approach.

Graphical interfaces are essential now for controling and monitoring the models of complex systems. On the 'output side' we have various tools, which use commercial programs to visualize the results. Currently, we are working in the direction of increasing the 'input side' to help the user define and set up the model system.

Appendix

To generate the random vector Dx_t the present GSA routine use an extension of the procedure used in reference [12]. We have calculated the Dx = g^-1(w) using a numerical integration of the visiting distribution probability gq_V(x). Where w is a random vector [0,1] obtained from an equi-probability distribution and g^-1 is the inverse of the integral of gq_V(x) given by

where q_A (q_V) is the acceptance index (visiting index), and

The integral of gq_V(x) has an analytical solution only if (q_A; q_V) = (1; 1) or (q_A; q_V) = (1; 2). For the general case it is necessary to make a numerical integration. In this paper, g^-1 has been calculated using a series integration and taking the inverse of a polynomial series, whose expansion is cut off at the 17^th order.

The integral in equation 17 can be written as

where

Solving equation 19 using power series, we have;

or

We look for an inverse function g^-1 = x = x(w) , such that x(w)-x = 0. From Eq.(21), we see that the inverse function may be expressed in terms of a power series,

This is guaranteed by a theorem: Let w(x) be analytic at x = x_o, and ¹ 0, then the inverse of w(x) exists and is analytic in a sufficiently small region about w(x_o) and its derivative is .

The coefficients B_n may be expressed in terms of A_n. However, it is possible to derive the coefficients more elegantly by employing Cauchy's formula. In this case, B_n takes the general form

The first few B_n coefficients can be expressed as

This last inversion was introduced in the GSA package, while the original proposition [11] uses a Lévy-Flight distribution as the visiting distribution.

References

[1]
• B.J. Alder and T.E. Wainwright, J. Chem. Phys. 27, 1208 (1957); Phys. Rev. 1, 18 (1970).
[2]
• N.L. Allinger and U. Burkert, "Molecular Mechanics'', ACS-Monograph, 177 (1982).
[3]
• E.P.G. Areas, P.G. Pascutti, S. Schreier, K.C. Mundim, P.M. Bisch, J. Phys. Chem. 99, 14882 (1995).
[4]
• E.P.G. Areas, P. G. Pascutti, S. Schreier, K.C. Mundim, P.M. Bisch, Brazilian J. Mol. Biol. Res. 27, 527 (1994).
[5]
• S. Kirkpatrick, J. Stat. Phys. 34 , 975 (1984).
[6]
• S. Kirkpatrick, C.D. Gelat and M.P. Vecchi, Science 220, 671 (1983).
[7]
• D. Ceperly, C. Alder, Science 231, 555 (1986).
[8]
• H. Szu and R. Hartley, Phys. Lett. A 122, 157 (1987).
[9]
• C. Tsallis, J. Stat. Phys. 52, 479 (1988).
[10]
• E.M.F. Curado and C. Tsallis, J. Phys A 24, L69 (1991); Corrigenda: J. Phys. A 24, 3187 (1991).
[11]
• C. Tsallis and D.A. Stariolo, Phys. A 233 , 395 (1996).
[12]
• K.C. Mundim, C. Tsallis, Int. J. Quantum Chem. 58, 373 (1996).
[13]
• M.G. Moret, P.G. Pascutti, P.M. Bish and K.C. Mundim, J. Comp. Chem. 19, 647 (1998).
[14]
• K.C. Mundim, T.J. Lemaire and A. Amin Bassrei, Physica A 252, 405 (1998).
[15]
• D.E. Ellis and J. Guo, in Electronic Density Functional Theory of Molecules, Clusters, and Solids, ed. D.E. Ellis, (Kluwer, Dordrecht, 1995) p263.
[16]
• L. Guo, D.E. Ellis, K.C. Mundim and B.M. Hoffman, Journal of Porphyrins and Phthalocyanines 2, 1, (1998)
[17]
• S. Dorfman and D. Fuks, Composite Interfaces 3, 431 (1996) and references therein.
[18]
• F. Dellanay, L. Froyen and A. Deruyttere, Journ. of Mater. Sci. 22, 1 (1987).
[19]
• U. Gangopadhyay and P. Wynblatt, Journ. of Mater. Sci. 30, 94 (1995); S. Hara, K. Nogi and K. Ogino, in Proceedings of Int. Conf. "High-temperature capillarity''.
[20]
• S. Dorfman and D. Fuks, Composites 27A, 697 (1996).
[21]
• D.E. Ellis, K.C. Mundim, D. Fuks, S. Dorfman and A. Berner, Mat. Res. Soc. Symp., Proc. 527, 69 (1998).
[22]
• H.W. Wyckoff, Acta Cryst. 11, 199 (1958).
[23]
• C-K. Guo, D.E. Ellis, O. Warschkow, unpublished.
[24]
• D.E. Ellis and K.C. Mundim, in Proc. of 9th CIMTEC-World Ceramics Congress Forum on New Materials, P. Vincenzine, Ed., Techna, Florence, 1998, to be published.
[25]
• P.G. Pascutti, K.C. Mundim, A. S. Ito and P.M. Bisch, J. Comp. Chem. (to be published 1999).
[26]
• M.A. Moret, M. Sc. Thesis (IF-UFBa 04/1996).

Publication Dates

History