California Institute of Technology
Title of the Effort
Highly Parallelized Large Scale Atomistic Simulations
for Design of Materials
Objective
The objective is to develop theoretical methodologies
for practical computations of the structures and
properties of real materials that can be used in industrial process
design for manufacturing new materials.
Principal Investigators
William A. Goddard III, wag@wag.caltech.edu
MMSC
Department of Chemistry, Beckman Institute (13974),
Caltech, Pasadena, CA 91125
Richard Freisner, rich@cucbs.chem.columbia.edu
CBS
Department of Chemistry, Columbia University, New York, NY 10027
ZhenGang Wang, zgw@macpost.caltech.edu
Department of Chemical Engineering (21041),
Caltech, Pasadena, CA 91125
Abhinandan Jain, jain@telerobotics.jpl.nasa.gov
Robotics, Jet Propulsion Laboratory (198219),
Pasadena, CA 91109
Stephen Taylor, steve@cs.caltech.edu
SCP
Department of Computer Science (25680),
Caltech, Pasadena, CA 91125
Related Information
The full 1995 report of the Materials and Molecular
Simulation Center in the Beckman Institute at Caltech is
also available on the World Wide Web (WWW).
This report
includes the detailed progress reports on the NSFGCAG project along with
projects supported by other agencies and by seven industrial
collaborators.
APPROACH
There is an enormous gap between current methodologies for
atomistic simulations and the level required for accurately describing
the relevant properties of industrially important materials. Our
strategy is to transcend from the most fundamental theory
(quantum mechanics, QM) to practical engineering designs in a sequence of
four or five levels as indicated in the following figure.
 Quantum Mechanics (QM)
First principles solution of the
Schrodinger equation (H F = E F) leads to accurate predictions of
properties (structures, chemical reactions, excited states); however, this
was limited to small systems (10 to 20 atoms) limiting the time and
distance scale to 1 ps and 1 nm.
 Force Fields (FF)
By averaging over the electrons from QM we
can obtain parameters (charges, force constants) while allowing
materials to be described in terms of atoms rather than
electrons.
 Molecular Dynamics (MD)
Using force fields (FF) the
fundamental equations become Newton's
equations (F = MA) rather than the Schrodinger equations. This allows
practical calculations on 1000 to 5000 atoms rather than
10 to 20, extending the time and distance scale to 1 ns and 10
nm.
 Coarse Graining (CG)
By averaging over the atoms for MD, we
can obtain parameters representing groups of
atoms (molecules, segments), considerably simplifying the calculations.
 Statistical Mechanics (SM)
Using the CG description we can
examine materials in terms of the large scale motions relevant to
macroscopic experiments. This extends the time and distance scales to
1 microseconds and microns.
 Continuum Parameters (CP)
The results from SM are combined
to obtain macroscopic or continuum
parameters (free energies, phase diagrams, partition
coefficients, solubility parameters) suitable for
practical chemical engineering software (e.g., ASPENPLUS) for design of
unit processes.
PROGRESS
Our team of 5 PI's has made significant progress on each of
the six tasks described under APPROACH. The new theoretical
methods developed for quantum mechanics and molecular dynamics that scale
sufficiently slowly with size so as to be practical for the very large
systems of industrial interest. Simulations of hundreds of atom systems for
quantum mechanics and tens of millions of atoms for complex systems
within the reach.
QUANTUM MECHANICS AND QUANTUM CHEMISTRY
 QM  FINITE.
For finite molecules; A new methodology (PSGVB) combining
pseudospectral (PS) multigrid and dealiasing strategies with
the sophisticated manybody wavefunctions [generalized
valence bond (GVB)] were implemented and applied to large
scale problems. PSGVB led to considerably better scaling with
size (N^2 rather than the N^3, N^4, N^5, N^6 characteristic of
alternative methods) and simpler parallelization.
parallelization. The basic PSGVB program has been
commercialized. PSGVB has been extended
to
 treat all atoms of the
periodic table (using core effective potentials),
 handle new
sophisticated wavefunctions (GVBRCI, MP2), and
 describe important properties (solvation
energies, hyperpolarizabilities).
Next we will (1) optimize these methods for
parallel implementations, (2) extend the
methodology to include GVBRCIMP2 and selfconsistent GVBRCI.
 QM  INFINITE.
Most practical materials
properties require a description of infinite systems using
periodic boundary conditions (PBC). This is threedimensional (3D) for bulk
properties or twodimensional (2D) for surface growth and interfaces. For
this purpose we have developed a new method, Gaussian dual space density
functional theory (GDSDFT), in which most parts scale linearly with N. In
implementing this we have developed a new
separable pseudopotential that can be applied to all
atoms of the periodic table.
We reformulated the theories for electronic structure calculations of
periodic systems in a way suitable for large scale calculations using
Gaussian basis functions. An accurate grid is introduced for efficient
calculation of matrix elements. A dualspace approach is used to calculate
the Coulomb potential with computational cost that scales linearly with the
size of basis set. A preconditioned generalized conjugate gradients
approach is introduced for rapidly converging wave functions expressed in
terms of Gaussian basis functions. This method is applied to a variety of
crystals (including diamond, GaN, AIN, CdTe, and C60) and surfaces
[including GaAs (110) and BN (110)] with excellent results.
The GDSDFT program is now being used for
several projects. The program optimized for Cray C90. For next stage we will
optimize it for parallel environments; KSR, SGI/PC, Intel/Paragon, Cray T3D and
IBM SP2.
FORCE FIELDS
The progress in force fields is twofold.
Over the past year new
specific force fields were developed utilizing ab initio quantum
chemistry programs with force field optimization programs.
These led to development of force fields for poly ethylene,
polysilane, nylons, polyoxymethylenes, polyvinyl (fluoride and
chloride) polymers, semiconductors and ceramics. For instance, the
Polysilane FF, denoted MSXX, was developed using the
Hessian biased method to describe accurately the vibrational states,
the ab initio torsional potential energy surface, and the ab
initio electrostatic charges of polysilane oligomers. This MSXX FF
was used to calculate various spectroscopic and mechanical properties
of the polysilane crystal. Stressstrain curves and surface energies
are reported. Gibbs molecular dynamics calculations (Nose,
RahmanParrinello) were used to predict various materials properties
at higher temperatures. Phonon dispersion curves and elastic constants
were calculated at various temperatures.
The Singular Value Decomposition method (SVD) has been implemented for
fast optimization of force fields using a variety of experimental and
theoretical data. In most instances, the experimental data is
incomplete or unavailable for a fully contrained force field
development. The SVD method provides an easy and automatic for force
field development for cases where parameters are not readily available
or where higher accuracy is necessary. Energy and Force evaluator
module of massively parallel program is enhanced to handle specific
forces necessitated by these force field developments.
MOLECULAR MECHANICS AND MOLECULAR DYNAMICS
The focus here is on extending the methods of MD to
physical systems of molecules, polymers, liquids, biopolymers and inorganic
materials with more than 10 million atoms per simulation cell while
accurately treating longrange interactions. The new
methodology is the cell multipole method (CMM). The CMM
methodology has been efficiently implemented into a general MD code on
the KSR and Intel/Paragon parallel supercomputers
(part of the Caltech NSFGCAG project). Both versions are in production
mode. Simulations carried out on both platforms for challanging materials
and biotechnology problems.
For fast internal coordinate dynamics on a million
atoms, we have developed the NewtonEuler Inverse Mass Operator
(NEIMO) method. This methodology has been extended to handle periodic
systems and has been incorporated into the advanced dynamics
(Gibbs/Nose) codes for single processor computers. Currently the technology
is being merged into MPSim Program environment using shared memory platform
as first test bed. The applications for this method will enable us to
investigate the long time dynamics of liquid polymer and solid interfaces;
which has tremendous impact on broad range of chemical industry applications:
wetting, adhesion, phase separation, coatings, dyes, etc.
Program is developed originally on KSR in an object oriented fashion to be
portable and efficient for other parallel computers. KSR version serves
as the basis for ports to parallel computers using a global memory model.
Developments using MPI on Cray T3D, IBM SP2 and SGI/PC is in the beta testing
stage.
Starting with the
th KSR code a communication layer was added for distributed memory
applications based on the active message strategy for the
JMachine. This has been implemented and tested on the Intel Touchstone
Delta/512 and Intel/Paragon at Caltech.
We now can do routine molecular dynamic calculations up to
1.5 million atoms with the KSR and Intel versions, and now are in production
for industrially important problems (diffusion, glass temperature)
for realistic description for amorphous polymers. A recent project is
launched in studying the aggregation of asphaltenes in crude oil which
presents great challenge to petrochemical industry using large scale
atomistic simulation of oil, asphaltene and water interface.
This constitutes a great leap
forward, allowing accurate atomistic simulations on systems 40 times
larger than previous methods.
Over the last few months, the molecular mechanics and
molecular dynamics program: Massively Parallel/Materials&Process
Simulation and Analysis (MPSim) Program System is
being find tuned for optimal performance
on various high performance computers:
 Shared memory MPSim on KSR and SGI Workstations;
 Message Passing MPSim on Intel using NX Express libraries;
 MPI MPSim on CRAY T3D, SGI Power Challenge, and IBM SP2
The program is the most general Molecular Mechanics and Dynamics program
available on massively parallel platforms:
 It can model any system: Periodic and nonperiodic structures, polymers,
peptides, nucleic acids, inorganic materials and organometallic
systems
 It uses generic force fields developed over the past ten years:
Dreiding II, Universal Force Fields, AMBER Force Fields, CHARMM Force Force
Fields and MM3 Force Fields.
 It can handle special force field developed using Biased Hessian and
FFOPT generated (Force fields generically known as MSXX)
 Long range interactions are uniquely and efficiently handled using
cell multipole methods (CMM and RCMM)
 A dynamic load balancing method is used to increase the efficiency.
 Robust Minimization Techniques are implemented
 Advanced Extended System Molecular Dynamics methods are included;
NoseHoover Canonical dynamics and ParinelloRahman Constant Pressure
dynamics.
 High frequency modes may be eliminated through rigid body dynamics with
Quaternion representation.
 Mixed mode dynamics for Polymers or biopolymers in solvents: rigid body
representation for solvents and stiff groups and fully flexible representation
for soft segments.
 NEIMO mechanics and dynamics: First consistent implementation of generalized
coordinate representation of many body newtonian mechanics and dynamics
methods developed for space science application to molecular mechanics and
dynamics of polymer and biopolymer simulations.
 Trajectory and structure analysis tools for simulations performed on
high performance computing platforms
 Visualization tools for trajectory and structure analysis.
 A graphical user interface running on desktop workstations
to prepare job control and job launching to high performance computers
 On line documentation of the program and utilities.
The first release/production version of the program is now available and
 used in studying
the thermodynamcis of
large scale atomic clusters which has profound effect in microcluster
physics.
 being used in studying large viruses: HRV14(0.5M atoms), and TBV
(3M atoms)
STATISTICAL MECHANICS
We have focussed on
developing systematic method for simulating mesoscopic and
macroscopic structural and thermodynamic properties of supramolecular
structures. This method enables us to efficiently sample phase space of
synthetic polymers and biopolymers to construct the partition function
from which the structural and thermodynamic properties may be obtained.
This method is being applied to dilute polymers melts, bulk polymer
systems, folding of medium size proteins.
Monte Carlo Program originally written in fortran is now being optimized
and rewritten in an object oriented fashion in C for compliance with
MPSim.
A new Monte Carlo method is developed and applied in estimating diffusion of
gases in glassy polymer networks. This
statistical method utilizes millions of step molecular dynamics (nanoseconds
long) run to explore a phenomena which is orders of magnitude larger than
the time scales accessible in Molecular Dynamics.
SIGNIFICANT EVENTS AND ACCOMPLISHMENTS
 PREDICTION OF GAS DIFFUSION IN POLYMERS.
A significant industrial application for atomistic simulations is the
prediction of diffusion of gases in glassy polymers. We are developing
a new approach where atomistic Molecular Dynamics is performed on
a large polymer system to extract compressibilities and the dynamical
behaviour of voids. This is then used in a Monte Carlo program to run
long trajectories for penetrant gas molecules. In the polymers studied
so far (PolyEthylene and PolyVinylChloride, PolyVinylyDeneChloride)
the qualitative trends and the quantitative diffusion constants agree
well with experimental data. This general method can be applied to any
polymer system with any penetrant diffusive molecule.
 PREDICTION OF DIELECTRIC PROPERTIES OF POLYMERS
Molecular dynamcis simulations at 300 and 400 K were used to calculate the
dielectric constants and dielectric loss for crystalline and amorphous
PVDF. These calculations used MSXX force fields developed from first
principles and validated with experimental structural and mechanical
properties of various forms of crystalline forms of PVDF. We showed that the
high dielectric constant for amorphous PVDF arises from the rapid changes in
the torsion angles during the Molecular Dynamics. These changes are
enhanced by solitonlike defects in chain that diffuse during simulation.
The dielectric losses lead initially to a stretched exponential decay
function.
 BAND OFFSETS OF HETEROJUNCTIONS
The use of localized Gaussian basis functions for large scale first
principles density functional calculations with periodic boundary
conditions (PBC) in 2 dimensions and 3 dimensions has been made possible by
using a dual space approach. This new method is applied to the study of
electronic properties of IIVI (II=Zn, Cd, Hg; VI=S, Se, Te, Po) and IIIV
(III=Al, Ga; V=As, N) semiconductors. Valence band offsets of
heterojunctions are calculated including both bulk contributions and
interfacial contributions. The results agree very well with available
experimental data. The p(2X1) cation terminated surface reconstructions of
CdTe and HgTe (100) are calculated using the local density approximation
(LDA) with twodimensional PBC and also using the ab initio
HartreeFock (HF) method with a finite cluster. The LDA and HF
results do not agree very well.
 LARGE HYPERPOLARIZABILITY PUSHPULL POLYMERS.
Recently significant advances have been made in engineering pushpull
organic chromophores to have very large hyperpolarizabilities (beta),
leading to materials with mu beta as high as 15000 X 10(48) esu. Such
developments have been slow and costy because of difficulties in synthesis,
purification, and measurement. As an alternative we have developed a new
quantum mechanical program (PSGVB/NLO) which provides predictions of beta
for such molecules far faster than previously possible. We have applied
PSGVB/NLO to predicting alpha, beta, and gamma for the high beta pushpull
organics and find excellent agreement with experiment. This suggests that
theory can be used as an effective tool for developing new nonlinear
optical materials.
 PREDICTION OF NEW DONORS FOR ORGANIC SUPERCONDUCTORS
The donors of all known one or twodimensional organic superconductors
are based on a core organic molecule that is either tetrathiafulvalene
(denoted as TTF) or tetraselenafulvalene (denoted as TSeF) or some mixture
of these two molecules. Coupling X, with appropriate accepters, Y, leads to
superconductivity. The oxidized form of X may be X(+) or X(2)(+) species in
the crystal. From ab initio quantum mechanical calculations
(HF/631G(**)),
we found that all known organic superconductors involve an X that deforms to
a boat structure while X(+) is planar. This leads to a coupling between
charge transfer and the boat deformation phonon modes that we believe is
responsible for the superconductivity of these materials. Based on this
idea we have developed similar organic donors having the same properties
and suggested that with appropriate electron accepters they will also lead to
superconductivity.
 MPSim MULTI PLATFORM PRODUCTION VERSIONS
In early fall, fully functional first production version of
Massively Parallel Materials and Process Simulation (MPSim) program is
completed. This version is a fully functional molecular mechanics and
molecular dynamics program for chemical, biologic and materials problems.
Its components are:
 User Interface: Tcl/Tk based graphical user interface
 Interface module: Input, output
 Force Field Evaluator module: Core of the program
 Simulators module:
 Energy  Force Evaluator
 Molecular Mechanics: Structure Optimizer
 Molecular Dynamics : NVE, NVT dynamics
 Analysis and Utility programs
 Documentation: Currently online, and tested by MSC members and associates.
 "Dielectric Properties of Poly(vinyledene fluoride) from Molecular
Dynamics Simulations,"
N. Karasawa and W. A. Goddard, III., Macromol. 6785 (1995).
 "Hessian Biased Force Field for Polysilane Polymers," C.B. Musgrave, S. Dasgupta and W.A.
Goddard, III, J. Phys. Chem. 99 13321 (1995).
 "Ab initio predictions of large hyperpolarizibility
pushpull polymers.
Julolidinylnisoxazolone and julolidinylnN,N'diethylthoabarbituric
acid", Daqi Lu, B.
Marten, Y. Cao, R. A. Freisner, and W.A. Goddard,
Chem. Phys. Lett. 242, 543, (1995).
 "Dualspace approach for densityfunctional calculations
of two and threedimensional crystals using Gaussian basis functions,"
X. Chen, JM. Langlois, and W.A. Goddard
Phys. Rev. B 52, 2348 (1995)
 "Effects of Catalyst Promoters on the Growth of SingleLayer
Carbon Nanotubes,"
CH. Kiang, W. A. Goddard III, R. Beyers, J. R. Salem, and
D. S. Bethune
Mat. Res. Soc. Symp. Proc. 359, 69 (1995)
 "First principles studies of band offsets at
heterojunctions and of surface reconstruction using Gaussian dualspace
density functional theory,"
X. Chen, A. Mintz, J. Hu,
X. Hua, and W.A. Goddard,
J. Vac. Sci. Technol. B13, 1715 (1995)
 "Building proteins from Calpha coordinates using the dihedral
probability grid MonteCarlo method," A. M. Mathiowetz and W. A. Goddard III,
Protein Sci. 4 , 12171232 (1995)
 "De novo prediction of polypeptide conformations using dihedral
probability grid MonteCarlo methodology,"
J. S. Evans, A. M. Mathiowetz, S. I. Chan, and W. A.
Goddard III
Protein Sci. 4 , 12031216 (1995)

"Stabilizing the boat conformation of cyclohexane rings,"
S. Dasgupta, Y.C. Tang, J.M. Moldowan, R.M.K. Carlson and W.A. Goddard,
J. Am. Chem. Soc. 117 , 65326534 (1995)
 "Prediction of new donors for organic superconductors,"
E. Demiralp and W.A. Goddard,
Synt. Metals 72, 297299 (1995)
 "Design and synthesis of a new peptide recognizing a specific 16basepair
site of DNA,"
C. Park, J.L. Campbell and W.A. Goddard,
J. Am. Chem. Soc. 117, 62876291 (1995)
 "Is carbon nitride harder than diamond  No, but its girth increases when
stretched (negative poisson ratio),"
Y.J. Guo and W.A. Goddard,
Chem. Phys. Lett. 237, 7276 (1995)
 " Experimental and theoretical studies of CO(CH4)(X)(+) with X=14,"
C.L. Haynes, P.B. Armentrout, J.K. Perry and W.A. Goddard,
J. Phys. Chem. 99, 63406346 (1995)
 "Free energy and surfacetension of arbitrarily large Mackay icosahedral
clusters,"
R.B. Mcclurg, R.C. Flagan, and W.A. Goddard,
J. Chem. Phys. 102, 33223330 (1995)
 "Parallel Calculation of ElectronTransfer and Resonance
Matrix Elements of HartreeFock and Generalized Valence Bond Wavefunctions,"
E. P. Bierwagen, T. R. Coley, and W. A. Goddard III
Am. Chem. Soc., 1994 ACS Symposium Series on Parallel
Computing in Computational Chemistry, T. G. Mattson editor, Chapter 7,
pg. 84 (1995).
FY95 PLANS
 PSGVB.
Parallelize the solvate,
hyperpolarizability, and GVBRCI methods.
Further develop and optimize the MP2, selfconsistent GVBRCI,
and GVBRCIMP2 methods. Apply to important industrial problems.
 GDSDFT.
Parallelize for the CRAY T3D. Apply to important industrial problems.
 Constant StressPressure Mechanics and Dynamics
Generalize the constant pressure and stress algorithms to handle
multimode mechanics and dynamics (for flexible, rigid and NEIMO threads)
Reformulate and fine tune the expensive stress calculations.
 Two dimensional periodicity for Interface Problems
Develop the interface simulation module to study extremely important
industrial processing and manufactruring problems revealing themselves
as interface issues: Wetting, friction, adhesion, wear, corrosion and
failure.
 Non Equilibrium Molecular Dynamics Module
Most of the materials problems reveal themselves as systems far from
equilibrium. Equilibrium molecular dynamics is insufficient to capture
the essense of the properties of the systems under the influence of
shear, or under thermal, pressure, concentration gradients; where system
shows a flow behavior. Synthetic hamiltonian formulations for performing
nonequilibrium molecular dynamics exist but applications have always
been restricted to simple model systems. The large scale massively
parallel platforms is making this step a feasible one. Especially,
the underlying general force field evaluator core of MPSim is our
strongest point to develop such a general NEMD program module.
TECHNOLOGY TRANSITION
Schrodinger Inc., Pasadena CA (contact Dr. Murco Ringnalda, ph:
8185689392, fax: 8185689778, email: info@psgvb.com)
has licensed the PSGVB technology from Caltech and Columbia. They
have a userfriendly interface, extensive user
guide, and extended the program in many ways. Over the past year they
made two commercial releases. PSGVB is also ported into Cerius2
interface (from Biosym/Molecular Simulations)
Chevron Petroleum Technology Co., La Habra, CA (Dr. Yongchun Tang,
ph: 3106947550) uses PSGVB, GDSDFT, and MDCMM to design new scale
inhibitors and corrosion inhibitors for oil
recovery applications (designing new and more effective calculations).
AlliedSignal, Morristown, NJ (Dr. Willis Hammond, ph:
2014554914 email: HammonB@mtomp201.research.allied.com) uses MDCMM
to predict structure and moduli of nylon polymers.
Asahi Chemical Corp. Japan (Terumasa Yamasaki email:
yamasaki@cs.fuji.asahikasei.co.jp, Takeshi Aoyagi, Shozo Yamanaka)
uses PSGVB, MDCMM to simulate a wide range of problems in catalysis,
electronic structure and material properties.
Asahi Glass Corp. Japan (Takashige Maekawa, Hiroshi Yamamoto, Takashi
Miyajima email: tam@agc.co.jp) uses MDCMM to simulate polmer properties.
Xerox Palo Alto Research Corp. (Dr. Ralph Merkle, email:
merkle@xerox.com) uses MDCMM to simulate nanosystems and understand
their vibrational properties.
Sandia National Labs (Drs. John Shelnutt and M. Cather Simpson email:
MCSIMPS@energylan.sandia.gov) uses MDCMM to obtain accurate force
fields for designing of porphyrin based chiral and achiral catalysts
for a variety of industrial problems.
Hughes Research Labs., Malibu, CA (Dr. Jenna Zinck, ph:
3103175913) uses GDSDFT to design metallorganics for MOMBE deposition of
HgCdTe
single crystal films.
Last modified on October 6, 1995.
Problems with this server should be reported to webmaster@www.wag.caltech.edu