Published: On

Molecular Dynamics Simulation of Crystalline C₆₀ Buckminsterfullerene with MedeA® and LAMMPS

This application note provides an overview of the forcefield based simulation of crystalline C₆₀ (Buckminsterfullerene) using the LAMMPS molecular dynamics simulation package. The emphasis is on the overall philosophy of LAMMPS calculations in the MedeA® environment.

LAMMPS and MedeA

LAMMPS1, the Sandia National Laboratories ‘Large-scale Atomic/Molecular Massively Parallel Simulator’ , provides impressive molecular dynamics performance, particularly when coupled with modern parallel-friendly compute environments.

LAMMPS has been employed in the simulation of systems comprising billions of atoms[^2] and may also be used in the detailed simulation of smaller systems.

Given LAMMPS wide adoption, and its parallel scaling attributes, there is keen interest in the convenient preparation of simulation models which may be employed with the program. This application note shows how LAMMPS and MedeA® work together.

LAMMPS commands allow for the construction of crystalline metallic systems, using a range of geometric manipulation commands. However, the creation of molecular systems is not specifically supported, and the creation of input datasets representing covalently bound systems for LAMMPS can be challenging, even for expert users, as the size, diversity, and complexity of molecular systems renders the model construction and in particular potential parameter assignment a time consuming and error prone undertaking.

Often hand written scripts are used to generateLAMMPS input data. However, script creation requires programming experience, and has a steep learning curve, which tends to reduce the diversity of molecular systems which are studied, as the effort required in creating input structures is significant.

MedeA® is the Materials Design® software environment for atomistic simulation and the MedeA LAMMPS interface provides for set up, calculation management, and analysis capabilities for LAMMPS users. Here we show an example of how a LAMMPS calculation can be conveniently prepared, administered, and analyzed using the Materials Design® MedeA® modeling environment.

## MedeA LAMMPS Interface

In order to perform molecular dynamics or molecular energy minimization calculations, LAMMPS requires a structural description of the system to be simulated, and details of the energy function and parameters (forcefield) to be employed in the calculation.

InfoMaticA returns more than 50 crystal structures for the element ‘C’, three of these structures are for the C₆₀ crystal.

Figure 1. InfoMaticA returns more than 50 crystal structures for the element ‘C’, three of these structures are for the C60 crystal.

Typically then, there are two key steps in the creation of structural input for LAMMPS: - model construction and - forcefield management.

## Model Construction Within MedeA® typical structural model construction methods include:

  • MedeA Crystal Builder.
  • MedeA PolymerBuilder.
  • MedeA Molecular Editor.
  • Supercell creation, Symmetry finding, and associated capabilities, as supported by the MedeA® modeling environment.
  • MedeA InfoMaticA. Provides access to hundreds of thousands of inorganic crystal structures, based on the ICSD, Pearson, Pauling, and NIST structural databases.
  • Specialized builder tools for the construction of Special Quasirandom Structures (SQS) alloy models, and nano-structures.
  • General crystalline model construction.
  • Editing operations (such as defect creation) applied to any of the models construction methods enumerated above.

    An illustration of the construction of a model based on the InfoMaticA database is provided in Figure 1. Here searching for the simple criterion of structures containing exclusively carbon within complete crystallographic information yields in excess of 50 structures. These represent diamond and graphite under different temperature and pressure conditions, and three of these hits are for structures of crystalline C₆₀. As an example, the structure reported by David and coworkers in 1991[^3] is selected here for further investigation.

     A single click takes you from the InfoMaticA database result to a three dimensional interactive model.

    Figure 2. A single click takes you from the InfoMaticA database result to a three dimensional interactive model.

    As Figure 2 shows, this crystal data readily yields its three dimensional structure which may be visualized in MedeA®.

    MedeA® supports a wide range of structural editing capabilities. In the present case, we are interested in simulating a larger system than that represented by the original crystal structure determination, so we create a P1, triclinic, supercell. To achieve this, we simply supply the number of unit cell repeats required in each crystallographic direction, and as shown in Figure 3, MedeA® constructs the appropriate supercell. MedeA® provides a wide range of surface, interface, and nano- builders, making structural model creation, even for complex systems, straightforward.

    Forcefield Management

    Once the structure has been built, LAMMPS simulations require suitable specification so that forcefield energy evaluation may be accomplished.

    This amounts to selecting a suitable set of forcefield parameters, identifying the atom types in the simulation structure and defining the functional forms and parameters to be employed in the LAMMPS calculation.

    MedeA® ships with a growing collection of forcefields, for systems ranging from inorganic glasses, to semiconductors, organic compounds, and metallic systems. In the present case, the PCFF+ forcefield has been employed. PCFF+ is based on the PCFF forcefield[^4] with extensions and refinements, particularly to non-bonded parameters, derived in a similar manner to those employed in the development of the COMPASS forcefield [^5] [^6] [^7] [^8].

    The association of forcefield parameters to particular sets of atoms within the structure is accomplished via a set of rules, which are contained in the forcefield file. Charges are assigned on the basis of bond increments which employ a representation of the polarization of chemical bonds to yield overall partial charges for the molecular systems described by the forcefield.

    Creating a 2x2x2 C60 simulation supercell

    Figure 3. Creating a 2×2×2 C₆₀ simulation supercell.

    Hence, based on atom type assignment rules and bond increment tables (both of which are supplied with MedeA) the C₆₀ structure is rapidly prepared for simulation with LAMMPS. Once the structure has been built, and the forcefield selected, just two mouse clicks are required to complete the atom type and partial charge assignment.

    Naturally, assigned forcefield types and energy functions may be inspected within MedeA®, in a convenient spreadsheet form. It is also straightforward to customize atom type and partial charge assignments, should the need arise, using the molecular spreadsheet.

    Constructing LAMMPS Simulations

    LAMMPS possesses a rich command language, which is thoroughly documented in the LAMMPS manual and on the LAMMPS website[^9].

    The [MedeA LAMMPS]/medea/LAMMPS) interface allows you to assemble visual flow charts describing complex simulation strategies in just a few operations. Figure 4 shows the essential approach to LAMMPS control employed by MedeA®. A set of simulation stages are linked to describe the overall simulation workflow.

    In the present case individual stages provide for the initiation of LAMMPS, the selection of non-bond interaction summation method, minimization of the cell, velocity assignment, and molecular dynamics property accumulation. The resulting flow chart may be shared, or applied to range of different simulation systems.

    Additionally, if custom LAMMPS commands, or scripting commands, are required to achieve a given modeling objective, these can be readily applied using the MedeA® interface, with appropriate custom stages in flowcharts. The MedeA LAMMPS interface makes routine calculations straightforward and yet provides for the full functionality of LAMMPS to be accessible in its interactive user interface.

    The interface allows you to construct LAMMPS simulations as flow charts

    Figure 4. The MedeA LAMMPS interface allows you to construct LAMMPS simulations as flow charts.

    Administering Calculations

    Monitoring and administering LAMMPS calculations is handled by the Materials Design® JobServer, as illustrated in Figure 5. Here you can check on the status of calculations, and if necessary, stop or restart calculations.

    The JobServer maintains all of the necessary information to reproduce calculation results, input structures, energy parameters, and LAMMPS input commands, structural, and energy function information, are stored with their associated LAMMPS output for future reference.

    This information is stored on disk in standard file formats, making it possible to easily document and reproduce calculations. Hence your time can be focused on applying simulation methods rather than bookkeeping calculation files.

    In the case of the present C₆₀ molecular dynamics calculation we can examine the input data files prepared by MedeA®. The structural and energy function data are shown in Figure 5. This file is constructed entirely automatically by MedeA® on the basis of the supplied structure and forcefield choice.

    The JobServer interface provides convenient access to all calculations performed on a given compute server.

    Figure 5. The JobServer interface provides convenient access to all calculations performed on a given compute server.

    The input files written by MedeA (see Figure 6, for example) are standard LAMMPS files and can be used with any desired compute server. Typically calculations are managed by MedeA’s JobServer, but, if required, the input files can be simply transferred to any desired compute server and calculations performed using the raw LAMMPS input information.

     Structural and energy function

    Figure 6. Structural and energy function written by MedeA® for use with LAMMPS.

    MedeA® writes the appropriate selection of energy functions for a given system to the LAMMPS input file. Each of the lines in the file shown describes individual energy terms for the C₆₀ structure, based on the PCFF+ parameterization.

    Analyzing Simulation Results

    Once a calculation has been completed, MedeA® provides a broad range of analysis capabilities. For example, as Figure 7 shows, the saved molecular dynamics trajectory may be animated and representative frames examined in detail. In the present case it is straightforward to create a movie showing the rotational disorder inherent in crystalline C₆₀, using video-capture software combined with MedeA’ s animation capabilities.

    Additional Capabilities

    MedeA® supports a range of forcefield types for LAMMPS simulations, including the class II, PCFF+ forcefield employed in the current example, Stillinger-Weber and Tersoff forcefields for semiconductor materials, and embedded atom method forcefields for the simulation of metallic systems. Additionally, MedeA/VASP provides for first-principles calculations, which may be employed in validating and complementing forcefield based simulations.

    MedeA Modules Employed in this Example

    The present calculations were performed with the MedeA® platform using the following integrated modules of the MedeA® software environment

  • The standard MedeA® framework including

    • crystal structure builders and
    • geometric analysis tools, as well as:
  • MedeA LAMMPS interface
  • MedeA JobServer
  • MedeA TaskServers
  • MedeA InfoMaticA with structural databases (ICSD, Pearson, Pauling, & NIST)

    Online Movie

    An animation showing 400 frames from the 100ps molecular dynamics calculation for crystalline C₆₀ at 200K is available at

 Interactively analyzing a LAMMPS trajectory for the C60 system.

Figure 7. Interactively analyzing a LAMMPS trajectory for the C₆₀ system.

Additional LAMMPS based molecular dynamics movies are may be found on the LAMMPS website:

### References

  1. S. Plimpton J. Comp. Phys. 117, 1 (1995) [^2]: M. Buehler, A. Hartmaier, M. Duchaineau, F. Abraham and H. Gao Acta Mech Sinica 21, 103 (2005) [^3]: R. Taylor, R. Ibberson, K. Prassides, D. Walton, J. Matthewman, J. Hare, W. David, H. Kroto and T.Dennis, Nature 353, 147 (1991) [^4]: H. Sun, S. J. Mumby, J.R. Maple, A. T. Hagler, J. Phys. Chem. 99, 5873 (1995) [^5]: H. Sun and D. Rigby, Spectrochimica Acta A153, 1301 (1997) [^6]: D. Rigby, H. Sun, B.E. Eichinger, Polymer International 44, 311 (1997) [^7]: H. Sun, J. Phys. Chem. B102, 7338 (1998) [^8]: H. Sun, P. Ren, J.R. Fried, Comput. Theor. Polymer Sci. 8, 229 (1998) [^9]: The LAMMPS website site