Structured programming is no longer enough for dealing with the large software projects allowed by today's computer hardware. An object-oriented computational model has been developed in order to achieve reuse, rapid prototyping and easy maintenance in large scale materials science calculations. The exclusive use of an object-oriented language is not mandatory for implementing the model. On the contrary, embedding Fortran code in an object-oriented language can be a very efficient way of fulfilling these goals without sacrificing the huge installed base of Fortran programs. Reuse can begin from one's old Fortran programs. These claims are substantiated with practical examples from a professional code for the study of the electronic properties of atomic clusters. Out of the about 20,000 lines of the original Fortran program, more than 70% of them could be reused in the C++ objects of the new version. Facilities for dealing with periodic systems and for scaling linearly with the number of atoms have been added without any change in the computational model.

We have investigated the macroscopic polarization and dynamical charges of some crystalline dielectrics presenting a mixed ionic/covalent character. First principles investigations have been done within the Hartree-Fock, LDA, and model GW approaches. All calculations have been performed on the same footing, using the all-electron FLAPW scheme. Apparently similar oxides have strikingly different behaviors: some (like the ferroelectric perovskites) have giant dynamical charges, while others (like ZnO) are quite normal and display dynamical charges close to the nominal static ones. We find the rationale for such differences.

Linear-response theory provides an efficient approach for calculating the vibrational properties of solids. Moreover, because the use of supercells is eliminated, points with little or no symmetry in the Brillouin zone can be handled. This allows accurate determinations of quantities such as real-space force constants and electron-phonon coupling parameters. We present highly converged calculations of the spectral function α2F(ω) and the average electron-phonon coupling for Al, Pb, and Li. We also present results for the free energy of vacancy formation in Al calculated within the harmonic approximation.

Determining the electronic and structural properties of semiconductor clusters is one of the outstanding problems in materials science. The existence of numerous structures with nearly identical energies makes it very difficult to determine a realistic ground state structure. Moreover, even if an effective procedure can be devised to predict the ground state structure, questions can arise about the relevancy of the structure at finite temperatures. Kinetic effects and non-equilibrium structures may dominate the structural configurations present in clusters created under laboratory conditions. We illustrate theoretical techniques for predicting the structure and electronic properties of small germanium clusters. Spefically, we illustate that the detailed agreement between theoretical and experimental features can be exploited to identify the relevant isomers present under experimental conditions.

We present an extension of a newly proposed tight-binding method of calculating total energies and related quantities to the study of binary compounds. As an example, we discuss the equation of state of the palladium-hydrogen system.

We present a tight-binding model which goes beyond the traditional two-center approximation and allows the hopping parameters and the repulsive energy to be dependent on the bonding environment. We show that this model works well for metallic as well as covalent systems.

We describe a new method that allows an efficient construction of high-quality initial wavefunctions which are required as input for iterative total-energy methods. The key element of the method is the reduction of the parameter space (number of wavefunctions) by about two orders of magnitude by projecting the plane-wave basis onto an atomic basis. We show that the wave functions constructed within this basis set are very close to the exact plane-wave wavefunctions, resulting in a rapid total-energy convergence.

Conjugate-gradient (CG) techniques for ab initio calculations of large complex systems have been examined for SIC-Al interfaces. The CG method by Bylander, Kleinman and Lee is more efficient than the Teter-Payne-Allan (TPA) method, the modified TPA method and the block Davidson method, although the TPA method is efficient for SiC surfaces. From the relaxed configurations, we have found strong attractive interactions between C and Al atoms, which should play a favorable role for adhesion between SiC and Al.

We have implemented a molecular-dynamics algorithm for silicon using a non-orthogonal tight-binding Hamiltonian with the functional form of Menon and Subbaswamy. Parameters for this Hamiltonian were determined by fitting to a database of first-principles total energy calculations of bulk phases and point defect formation energies. These geometries were chosen to reproduce the configurations seen in defective crystalline and amorphous silicon. We have also implemented the non-orthogonal density-matrix method, paying particular attention to data motion locality to facilitate efficient parallelization of the algorithm. The necessary sparse matrix operations (trace, transpose, matrix multiplication) have also been implemented on a single processor workstation with an algorithm which takes O(N) time. Tests of the method's accuracy involved calculations of surface energies and structural reconstructions and activation energies for bulk diffusion through concerted exchange. We present results of a simulation of the melting and rapid quenching of a silicon sample using molecular-dynamics, and examine the resulting structures.

A semi-empirical methodology [1,2] developed to model and simulate covalently bonded networked systems is modified to study the heteroatomic mixtures of silica. This methodology is capable of grasping the essential qualitative and quantitative features of the coupling between the electronic coordinates and the geometric structure. The methodology is used to simulate and to probe the structural and thermodynamic properties of the bulk crystalline, amorphous solid and the melt states of silica.

In a previous paper [1], we developed a form of variational trial wave function and applied it to van der Waals clusters: five or less atoms of Ar and Ne modeled by the Lennard-Jones potential. In addition, we tested the trial functions for a hypothetical, light atom resembling Ne but with only half its mass. We did not study atoms such as He4 with larger de Boer parameters, i.e., systems in which the zero point energy plays a more important role relative to the potential energy. This is the main purpose of the present paper. In fact, we study clusters to the very limit where the zeropoint energy destroys the ground state as a bound state. A simple picture of this un-binding transition predicts the power law with which the energy vanishes as the de Boer parameter approaches its critical value and the power of the divergence of the the size of the clusters in this limit. Our numerical results are in agreement with these predictions.

We have implemented the fully relativistic and spin-polarized extension of the O(N) real space locally self-consistent multiple scattering formalism [1]. Here we present preliminary results for the spin and orbital magnetic moments and magneto-crystalline anisotropy in iron metal. The spin and orbital moments are in good agreement with the results of conventional electronic structure methods. In addition we present preliminary results for the magnetocrystalline anisotropy obtained in the frozen potential approximation.

An analytic procedure for inverting ab initio cohesive energy vs. volume curves to obtain a parameter-free interatomic potential is presented. The procedure determines the radial functions in a cluster potential, assuming an angular dependence. The method is a nonlinear generalization of the ab initio pair potential formula of Carlsson, Gelatt and Ehrenreich to higher orders of cluster expansion [1]. We demonstrate our method by deriving an ab initio cluster potential for silicon with the Stillinger-Weber angular dependence [2]. A novel property of the potential is that bond-bending forces are not monotonic in the bond length, as is the case with most empirical potentials, but rather are peaked at the first neighbor distance. The validity of the inverted potential is discussed along with its use in the development of new interatomic potentials.

We present a method for selfconsistent Density Functional Theory calculations in which the effort required is proportional to the size of the system, thus allowing the aplication to problems with a very large size. The method is based on the LCAO approximation, and uses a mixed approach to obtain the Hamiltonian integrals between atomic orbitals with Order-N effort. We show the performance and the convergence properties of the method in several silicon and carbon systems, and in a DNA periodic chain.

We briefly review recently developed fast algorithms for the evaluation of electrostatic fields in composite materials consisting of a collection of piecewise homogeneous inclusions embedded in a uniform background. These algorithms are based on combining a suitable boundary integral equation with the fast multipole method and a conjugate gradient-like iterative method. The CPU time required grows linearly with the number of points in the discretization of the interface between the inclusions and the background material, bringing large-scale calculations within practical reach.

In our continuing efforts towards designing materials with controlled optical properties, largescale molecular dynamics simulations of a molecular cluster of a liquid crystalline cyclic siloxane are still limited by the size of the molecular system. Such simulations enable evaluation of the orientation order parameter of the system, as well as modelling the behavior of the material in bulk. This study summarizes improvements in the implementation of the fast multipole algorithm for computing electrostatic interactions which is included in the molecular dynamics program PMD[7, 8], such as the elimination of computations for empty cells and the use of optimal interaction lists. Moreover, an improved implementation of a 3-D Fast Multipole Method (FMM3D) based on the algorithm previously proposed[1, 2] is described in detail. The structure of the module, details of the expansions, parallelization, and its integration with the molecular dynamics simulation code are explained in detail. Finally, the utility of this approach in the study of liquid crystalline materials is briefly illustrated.

We introduce an efficient and scalable parallel implementation of tight-binding molecular dynamics (TBMD) which employs reordering of the atoms in order to maximize datalocality of the distributed tight-binding (TB) Hamiltonian matrix. Reordering of the atom labels allows our new algorithm to scale well on parallel machines since most of the TB hopping integrals for a given atom are local to the processing element (PE) therefore minimizing communication. The sparse storage format and the distribution of the required eigenvectors reduces memory requirements per PE. The sparse storage format and a stabilized parallel Lanczos eigen-solver allow consideration of large problem sizes relevant to materials science. In addition, the implementation allows the calculation of the full spectrum of individual eigen-values/-vectors of the TB matrix at each time-step. This feature is a key issue when the dielectric and optical response must be computed during a TBMD simulation. We present a benchmark of our code and an analysis of the overall efficiency.

Using a reactive empirical bond-order potential (REBOP) model for hydrocarbons1, large scale molecular dynamics simulations of carbon systems are carried out on parallel machines. Structural and dynamical correlations of amorphous carbon at various densities are studied. The calculated structure factor agrees well with neutron scattering experiments and the results of tightbinding molecular dynamics simulations. The dynamic behavior of crack propagation through graphite sheet is also investigated with the molecular-dynamics method. Effects of external stress and initial notch shape on crack propagation in graphite are studied. It is found that graphite sheet fractures in a cleavage-like or branching manners depending on the orientations of the graphite sheet with respect to the external stress. The roughness of crack surfaces is analyzed. Two roughness exponents are observed in two different regions.

We discuss the implementation of classical molecular dynamics simulations on Single- Instruction Multiple-Data computers with two-dimensional mesh connectivity. Specific timing results are given for the MasPar MP-1, complemented by a general time complexity analysis. Our main interest is in simulations with a number of particles N that is comparable to the number of processors P. In particular, we discuss results for large Si-clusters with up to 2000 particles, using a Stillinger-Weber potential. Particular attention is paid to the problem of handling in an efficient way the calculation of the three-body potential and force.

We present results of ab-initio electronic structure calculations and molecular dynamics simulations of organic molecules carried out using adaptive curvilinear coordinates, within the local density approximation of density functional theory. This approach allows for an accurate treatment of first-row elements, which makes it particularly suitable for investigations of organic compounds. A recent formulation of this method relies on a real-space approach which combines the advantages of finite-difference methods with the accuracy of adaptive coordinates, and is well suited for implementation on massively parallel computers. We used molecular dynamics simulations to obtain the fully relaxed structures of nitrosyl fluoride (FNO), and of the aromatic heterocycles furan and pyrrole. The equilibrium geometries obtained show excellent agreement with experimental data. The harmonic vibrational frequencies of furan and pyrrole were calculated by diagonalization of their dynamical matrix and are found to agree with experimental data within an rms error of 25 cm-1 and 28 cm-1 for furan and pyrrole respectively. This accuracy is comparable to that attained for smaller organic molecules using elaborate quantum chemistry methods.