SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
Integrating Open-Source Software Applications to Build Molecular Dynamics Systems Bruce M. Allen,* Paul K. Predecki, and Maciej Kumosa Three open-source applications, NanoEngineer-1, packmol, and mis2lmp are integrated using an open-source file format to quickly create molecular dynamics (MD) cells for simulation. The three software applications collectively make up the open-source software (OSS) suite known as MD Studio (MDS). The software is validated through software engineering practices and is verified through simulation of the diglycidyl ether of bisphenol-a and isophorone diamine (DGEBA/IPD) system. Multiple simulations are run using the MDS software to create
MD cells, and the data generated are used to calculate density, bulk modulus, and glass transition temperature of the DGEBA/ IPD system. Simulation results compare well with published experimental and numerical results. The MDS software prototype confirms that OSS applications can be analyzed against real-world research requirements and integrated to create a C 2014 Wiley Periodicals, Inc. new capability. V
Introduction
calculating physical properties of an epoxy material. This was achieved by calculating the density of the diglycidyl ether of bisphenol-a and isophorone diamine (DGEBA/IPD) epoxy system over a large temperature range allowing the simulated epoxy to transition from the glassy to viscous state. The glass transition temperature and the bulk modulus of the simulated epoxy were also calculated. The DGEBA/IPD system was selected for simulation due to its complexity and the fact that this system had been previously investigated providing experimental and numerical data for comparison.[7] The DGEBA and IPD structures are shown in Figures 1 and 2, respectively. The actual epoxy material consists of a three-dimensional network of polymer and cure molecules intertwined with other networks of itself creating a crosslinked and intertwined amorphous system. For simplicity, oligomers composed of nine DGEBA and four IPD molecules were used in the simulation in an attempt to model the epoxy material. Each oligomer was 75% intracrosslinked within itself as shown in Figure 3. A total of five oligomers were used to create a cell, which provided intertwining to the cell. Tack and Ford modeled and simulated a similar system consisting of DGEBF and DETDA using the Amorphous Cell in Materials Studio by Accelrys.[8,9] Modeling a similar system using MDS was important due to the intertwining of oligomers provided by packmol.[6] The Consistent Family of Force Fields (CFF91) was used for all simulations.[10] Figure 4 illustrates through a flowchart the three components of the MDS software applications, shown in yellow, and the LAMMPS application shown in green. The integration between the three applications was accomplished with
The goal of this work was to integrate existing open-source software (OSS) applications for the purpose of quickly creating an initial molecular dynamics (MD) cell, also known as a system of molecules or a system for simulation using the largescale atomic/molecular massively parallel simulator (LAMMPS).[1] LAMMPS and the visual MD (VMD) software are used to run MD simulations and to view atom trajectories.[2,3] LAMMPS itself is an object toolkit and can be built as a library and called as part of a larger software application. Colleagues reported initial epoxy system creation times of up to 6 months due to the complexity of specifying single system geometry and associated force field data for LAMMPS. Also, ring spearing and manually fixing systems generated with commercial software made finding alternative system creation approaches attractive.[4] A search for chemistry or nanotechnology computer-aided design (CAD) software yielded NanoEngineer-1 (NE-1), while a review of tools bundled with LAMMPS provided msi2lmp.[5] The initial project direction was to integrate the two software applications, but NE-1 lacked a MD cell or cell concept and drawing medium to large geometries (at least several thousand atoms) would take significant time. Later, while attending the LAMMPS Workshop in August 2011, packmol was brought to our attention, which had the necessary cell concept explicitly designed into it.[6] Packmol was determined to be the missing software application needed to complete the OSS development. The three integrated applications, NE-1, packmol, and msi2lmp are now known collectively as MD Studio (MDS). Initial development requirements were based on the need for application source code and a development environment to support software modifications and test. A lengthy comparison of many OSS applications was not attempted.
Simulated Epoxy System The MDS software integration was tested in this research by following standard software engineering practices and by 756
Journal of Computational Chemistry 2014, 35, 756–764
DOI: 10.1002/jcc.23537
B. M. Allen, P. K. Predecki, M. Kumosa University of Denver, Daniel Felix Richie School of Engineering & Computer Science, Mechanical & Materials Engineering, 2390 South York Street, Denver, Colorado 80210 E-mail: [email protected] Contract grant sponsor: National Science Foundation [(NSF) Grant Opportunities for Academic Liaison with Industry (GOALI)]; contract grant number: #CMMI-1232520 C 2014 Wiley Periodicals, Inc. V
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
Figure 1. Diglycidyl ether of bisphenol a (DGEBA) epoxy monomer resin.
Figure 2. Isophorone diamine (IPD) crosslinker (curing agent).
the development of the enhanced molecular machine part (EMMP) file. The EMMP file was developed for sharing data among the three applications, and is based on the existing molecular machine part (MMP) file. The LAMMPS geometry input file (LGIF) is the final output from the MDS software. The LGIF and the LAMMPS Command Input File were provided to the LAMMPS software, and the LAMMPS output trajectory file was generated and used to calculate the final simulation results (e.g., Tg, density, volume, temperature) for the epoxy system of interest. The trajectory files generated by LAMMPS were viewed using the VMD application.
Software Use Following the flowchart in Figure 4, it took approximately 2 h to draw a single oligomer. Single DGEBA and IPD molecules were drawn and then copied several times. Subsequently, the molecules were attached to one another to form the oligomer.
Figure 3. Oligomer structure.
SOFTWARE NEWS AND UPDATES
Care was taken to not introduce high energies into the starting molecules or when attaching them to form the oligomer. The NE-1 molecular mechanics energy minimzer (MMEM) was executed periodically to lower the molecules’ and oligomers’ potential energies. The existing NE-1 “combine” feature was used to group the molecules in the oligomer into a single molecule. Grouping the attached molecules into a single molecule (oligomer) was necessary because the EMMP file parser added to packmol, expected a single molecule per file. The final EMMP file created was a template for packmol containing the oligomer. The initial cell of five oligomers is shown in Figure 5. Running packmol is an iterative process; therefore, an initial guess of the initial cell size and shape was made based on the size of the oligomer and the number of oligomers. Packmol was run several times, and the volume of the cell was reduced until no solution was obtained. This previous successful execution provided the minimum volume for the initial cell. Packmol took no more than 1 h to generate the initial cell. Finally, the LGIF was created by msi2lmp in a few seconds. The DGEBA/ IPD system required no more than 3 or 4 h to generate the initial cell. The CFF91 force field file required maintenance to add generic values for missing parameters reported by msi2lmp. Maintenance times are outside the execution of the MDS software, but do have an obvious affect on the simulation results.
Software Features and Modifications NE-1 is CAD software created by Nanorex for molecular-level CAD.[5] The source code consists of Python and C. The WingIDE integrated development environment and debugger was used to test Python software changes. NE-1 is capable of modeling and simulating nanomachines, DNA, carbon nanotubes, and many other systems both organic and inorganic. NE-1 stores complex geometry in its MMP file format. The MDS software uses the MMP file as a molecule template. It provides the positions of atoms and bonding topology from which bond lengths, angles, and dihedral angles can be calculated. The MMP file format was expanded in this research to include atom types, thus creating the EMMP file, and was used for building the DGEBA/IPD initial cell. The original NE-1 uses a force field, NanoDynamics-1 (ND-1), which is element based and cannot be used for MD simulations in LAMMPS. ND-1 provides a single atom type for all chemical environments, whereas most force fields used with LAMMPS provide several per atom. Therefore, the modifications to the MMP file were critical. Journal of Computational Chemistry 2014, 35, 756–764
757
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
Figure 4. Flowchart for MDS software.
Modifications to the NE-1 graphical user interface (GUI) were made to support user access and selection of atom types while drawing molecules. Figure 6 shows the force field chooser on the left hand side of the NE-1 GUI. The atom chooser shows oxygen highlighted in blue, and the generic
oxygen atom type is highlighted in the force field chooser. Notice the CFF91 force field is selected in Figure 6. On the drawing canvas, an oxygen atom is highlighted in yellow, to allow the text and numeric atom types to be displayed. In the interest of completing the prototype quickly, atom types are given numerical representations to preclude token parsing issues and allow every element to have a maximum of 100 atom types. The ToolTips dialog is displayed in Figure 7. The dialog allows the user to select whether numerical or character atom types are displayed when atoms are highlighted on the canvas as presented in Figure 6. The user can choose another force field from the drop down box in the force field chooser. Clicking on the CFF91 title button displays another force field to be selected. The COMPASS atom types, not the associated proprietary force field parameters, were integrated into the NE1 force field chooser as a proof of concept, demonstrating how multiple force fields can be integrated into the application.[11] Packmol was developed by Dr. Leandro Martinez at IMECCUNICAMP.[6] Packmol was written in FORTRAN and the GNU FORTRAN compiler and the dynamic data display debugger were used in the CYGWIN environment under Windows 7 to make and test source code changes. The packmol application reads one or more molecule templates and fills a user defined volume with a user specified number of each molecule. Packmol allows the user to generate many unique cells by specifying a unique pseudorandom number seed. This allows the user to create and later recreate a cell based on initial system configuration information.[6] The original packmol software could not read or write MMP files, and was, therefore, modified to read and write the EMMP file. An example of how packmol functions is as follows: a cube containing 1,000 DGEBA molecules and 500 IPD molecules is created using a single-molecule template for DGEBA and single-molecule template for IPD. Given the size of the DGEBA and IPD molecules and the packing constraints, intertwining of molecules is not expected. This is not the case for the oligomers used in this research.
Figure 5. Initial cell: Five oligomers. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
758
Journal of Computational Chemistry 2014, 35, 756–764
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
SOFTWARE NEWS AND UPDATES
Figure 6. New software feature allows choosing the force field. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Figure 7. New software feature allows display of numeric atom types. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
The oligomers used in this work were unexpectedly intertwined, which was an added benefit to cell creation. For the cell shown in Figure 5, packmol calculated a cubedshaped initial cell with ˚. a side length of 55 A The initial cell was then reduced to the equilibrated final cubedshaped geometry, 30.5 ˚ per side at 1 atm and A 298 K using LAMMPS. The final cell in Figure 8 is rendered in unwrapped atom coordinates.[2] The simulation was run using periodic boundary conditions (PBC) which wrap atoms that leave the boundaries of the cell into positions directly on the opposite side of the cell. VMD was used to display the final cell. VMD attempts to display elements using unique colors. Figure 8 illustrates the point that the final cell is considerably denser than the initial cell seen in Figure 5. Packmol optimized the packing of the molecules by calculating a minimum cost function value. The minimum cost is associated with molecules being very close to one another and providing a minimum system volume, while minimizing the repulsion forces between atoms of different molecules. The cost function contained a user defined constraint for maintaining
Journal of Computational Chemistry 2014, 35, 756–764
759
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
Table 1. MDS requirements map. Requirement Experimental density Crosslinked and intertwined cell geometry/moduli Minimize repulsive forces Force field support Calculate atom partial charges Draw molecules Manually type atoms Modify and add atom types Generate LGIF Uniquely identify molecules
NE-1
packmol
Msi2lmp
LAMMPS
X X
X X
X
X
X
X X
X X
X X
X X X X X
The density requirement determines the free volume in the epoxy cell. Crosslinking and intertwining of polymer chains determine moduli of the material and density. The remaining requirements were cited above as software modifications to the three software applications. Table 2 lists the specific source code changes by application and programming language. Figure 8. Final cell: five oligomers after simulated annealing. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
a minimum distance between atoms of different molecules ˚ . When atoms from different molewith an initial value of 2 A cules were too close to one another, high repulsive forces occurred, resulting in atoms moving a large distance in a single simulation time step. LAMMPS terminates under such conditions due to bad dynamics. LAMMPS supports identification of individual molecules to calculate mean-squared displacement (MSD) during simulation. Packmol was modified again to uniquely identify molecules in the EMMP file allowing MSD calculations for small molecules in the cell. The molecule ID number generated in packmol was passed to msi2lmp for inclusion in the LGIF. The msi2lmp software was created by Dr. John Carpenter at Cray Research.[2] Msi2lmp generated a LGIF from the Cartesian coordinate file/molecular data file (CAR/MDF) file pair, which is one file format supported by Materials Studio, a commercial software application.[9] Msi2lmp was written in the C language, and the Microsoft Visual Studio tool was used to modify and test its source code modifications. The msi2lmp software enumerates all unique bond and nonbonded force relationships among atoms of a molecule and between many molecules composing a system. Msi2lmp also populates the force field parameters from a user supplied force field file into the LGIF. The epoxy material simulation required atom partial charges for Coulombic force calculations in LAMMPS. The original msi2lmp transferred the atom partial charges from the (CAR/MDF) file pair directly to the LGIF when processing files created by Materials Studio. The modified msi2lmp software, in this research, calculates the atom partial charges using the method defined by Sun after reading atom data from the EMMP file.[11] The LAMMPS “full atom” model is written to the LGIF, which is required to support the atom partial charges and the molecule ID created by packmol.[2] The above software modifications reflect the new requirements needed to integrate the three applications. Table 1 summarizes the MDS requirements. 760
Journal of Computational Chemistry 2014, 35, 756–764
Software Verification The MDS software developed in this study was subsequently verified by using standard software test methods of white and black box testing, requiring the use of debuggers to examine data values and pathways taken in the modified software. The MDS software was parallel tested using a single DGEBA/IPD oligomer to generate two LGIF files using two separate paths through the modified msi2lmp software. Parallel software testing was performed by converting the EMMP file into a protein data bank (PDB) file, which was then transformed into a CAR/MDF file pair. The CAR/MDF file part was subsequently converted to a LGIF using the existing pathway in the msi2lmp application for string-based atom types. The EMMP file was rewritten into a LGIF using a different pathway in the msi2lmp application that understood numeric atom types. The Beyond Compare software tool was used to compare the two LGIF files for differences.[12] Figure 9 shows the process for the above parallel testing. The two LGIF files were nearly identical, thus showing that the new MDS software created the correct file for simulation. Subtle differences in the two LGIF files were traced back to how NE-1 created errors in the PDB file. The MD verification simulations were executed using the major simulation parameters listed in Table 3.
Table 2. Software modifications. Software modifications
NE-1
NE-1 EMMP file parser (Python) Packmol EMMP file parser (FORTRAN) Msi2lmp EMMP file parser(C) GUI added force field chooser (Python) GUI changes to tooltips dialog (Python) Passes atom types Calculate atom partial charges (C)
packmol
Msi2lmp
X X X X X X
WWW.CHEMISTRYVIEWS.COM
X X
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
Table 4. Single round of simulated annealing at 1 atm. Step 1 2 3 4
Start temperature (K)
End temperature(K)
Duration(ps)
298 298 600 600
298 600 600 298
50 50 200 50
Begin temperature (K)
End temperature(K)
Duration (ps)
298.0 298.0 328.2 328.2 358.4 358.4 388.8 388.8 419.0 419.0 449.2 449.2 479.4 479.4 509.6 509.6 539.6 539.6 570.0 570.0 600.2
298.0 328.2 328.2 358.4 358.4 388.8 388.8 419.0 419.0 449.2 449.2 479.4 479.4 509.6 509.6 539.6 539.6 570.0 570.0 600.2 600.2
10 100 10 100 10 100 10 100 10 100 10 100 10 100 10 100 10 100 10 100 10
Table 5. LAMMPS NPT ensemble. Step
Figure 9. Parallel testing msi2lmp. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
PBC were used to maintain a constant number of atoms over the simulation duration and to isolate the computations from surface effects. The Van der Waals and Coulombic cutoffs ˚ to support the LAMMPS full atom were set to 9.5 and 10.0 A model. The commonly used in organic system simulations, 1 fs time-step, was assumed.[13]
Simulation Methods and Results The following simulation steps were performed to generate the final cell: MMEM, quench (velocity rescaling), and simulated annealing (SA). Table 4 presents the SA details for each of the five rounds. System equilibrium, after SA, was accomplished by running a series of NPT ensembles, where NPT represents holding the number of moles (N), pressure (P) and temperature for the system constant while allowing the system volume to change over time. The ensembles are listed in Table 5. The bulk modulus was calculated by linearly increasing the system pressure from 1 to 5001 atm at 298 K using the NPT ensemble in Table 6. Table 7 lists the calculated density values based on cell dimensions and temperature. The cell dimensions were
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
collected every 0.1 ps, and the calculated density at that temperature was the average of 1,000 data points. The simulated density value was verified by comparing it to the experimentally determined value from the literature for the same resin system tested at ambient conditions (i.e., 298 K and 1 atm).[10] The agreement was within 7% as shown in Table 8. This clearly indicates that MDS provided a high-quality estimate of the resin’s density. The cell dimensions and the bulk modulus of the resin at ambient temperature are presented in Tables 9 and 10, respectively. In addition, a comparison is made in Table 10 between the experimental and simulated values of the modulus. Similar to the density predictions, the error in the modulus prediction is small and was found to be less than 3%. MDS was further verified by determining the glass transition temperature of the resin of interest. Figure 10 presents a graph of the average density data plotted, in black diamonds, against the average temperatures from Table 7. The line in black is the least squares fit of the first five temperatures vs. density values.
Table 3. LAMMPS simulation parameters. Parameter PBC K-space Solver Van der Waals cuffoff Coulombic force cutoff Time step
Value Default PPPM 9.5 A˚ 10.0 A˚ 1 fs
Table 6. Bulk modulus data collection at 298 K NPT ensemble. Step 1 2 3
Begin pressure (atm)
End pressure (atm)
Duration (ps)
1.0 1.0 5001.0
1.0 5001.0 5001.0
100 100 100
Journal of Computational Chemistry 2014, 35, 756–764
761
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
Table 7. Simulated densities by NPT ensemble. Temperature (K)
Lx (A˚)
˚) Ly (A
˚) Lz (A
Volume (cc) 3 1.0e20
Density (gm/cc)
297.55 328.02 358.12 388.63 418.66 449.02 479.55 509.85 539.95 570.57 599.86
30.75 30.74 30.78 30.83 30.79 30.93 31.07 31.39 31.32 31.46 31.56
30.76 30.76 30.79 30.84 30.80 30.94 31.08 31.14 31.34 31.47 31.57
30.77 30.77 30.80 30.85 30.81 30.95 31.09 30.42 31.35 31.48 31.58
2.8954 2.9106 2.9197 2.9344 2.9222 2.9634 3.0042 3.023 3.0783 3.1164 3.14645
1.056 1.0510 1.0458 1.0427 1.0470 1.0325 1.0185 1.0122 0.99404 0.98187 0.9725
The remaining points were used to generate the least squares fit of the remaining six temperatures vs. density values as indicated by the gray dashed line. The first five points were selected to represent the glassy material prior to Tg, whereas the remaining points represent the viscous material. Tg was also estimated by selecting the first six and four glassy points as shown in Figures 11 and 12. Taking the average value of three determined Tgs using the first four, five, and six glassy data points, the average Tg was predicted to be 433 K.
Discussion Particle–particle–particle-mesh (PPPM) performs nonbonded calculations faster than using the Ewald method in a multicore environment. It fails sooner in a multicore environment, due to missing atoms compared to the Ewald method, making PPPM preferable.[2] The results shown in Table 7 took approximately 3 days on a 2.93 GHz 4-core i7 Intel processor using the PPPM option, whereas the same simulation required approximately 5 days using the Ewald option. Using a combination of MMEM, quench, and SA allowed cell conformations to reach a density near the experimental value at ambient conditions. It took approximately 36 h to achieve the ambient density shown in Table 7. Simulation of a system of small molecules would not have required such a long equilibration time. The oligomer size and perhaps the intertwining are responsible for the long preparation time. The oligomers were used to achieve intracrosslinking in the epoxy, due to previous attempts to manually intercrosslink oligomers led to longer equilibration times. The longer times were due to the introduc-
tion of nonequilibrium geometry, when constructing the crosslinks using NE-1. The calculated system density undershoots by 7% compared to the experimental value.[10] Polymer density is also a function of the degree of crosslinking.[14] Shokuhfar and Arab calculated DGEBA/DETA densities for crosslinking densities from 5 to 81%. The experimental values for 16, 25, and 37% are 1.16, 1.19, and 1.13 gm/cc, respectively, whereas the simulated densities are 1.1, 1.11, and 1.12 gm/cc, respectively.[15] In general, simulations undershoot the experimental values. Epoxy materials absorb water vapor, increasing density. There appears to be enough uncertainty in the experimental and simulation methods making determination of monotonically increasing densities with increasing degree of crosslinking difficult. The fifth density in Figure 10 is higher than expected and was most likely due to the heating rate and equilibration time used in elevating the system temperature from 388 to 418 K. The system is transitioning from glassy and viscous, and further
Table 10. Simulated bulk modulus.
D (MPa)
DV (cc)
V0 (cc)
Bulk modulus (GPa)
Experiential bulk modulus (GPa)
506.6
2.90e-21
2.94e-20
5.13
5.0110
Table 8. Simulated vs. experimental density at 1 atm and 298 K. Simulated density (gm/cc)
Experimental density (gm/cc)
Percent difference
1.056 6 0.02
1.13110
761
Table 9. Simulated bulk modulus raw data.
762
Pressure (atm)
˚) Lx (A
Ly (A˚)
˚) Lz (A
Volume (cc)
1.0 5001.0
30.86 24.76
30.87 29.82
30.88 35.92
2.94e-20 2.65e-20
Journal of Computational Chemistry 2014, 35, 756–764
Figure 10. Five point glassy Tg least squares.
WWW.CHEMISTRYVIEWS.COM
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
moleculardynami/, and is chronicled on our blog at http:// moleculardynamicsstudio.blogspot.com/. Several software enhancements are in progress including static and dynamic crosslinking to build better cells, automatic atom typing, support for all open-source atomistic and united atom force fields, and integration of packmol and msi2lmp applications into NE-1 for ease of use. Better crosslinking will result in ease of use and speed up preparing the final cell for simulation. Automatic atom typing will free the user to manually type exceptional cases, where the software is unable to determine the correct atom type. The CFF force fields usually include an “auto equivalence” table to provide missing force field parameters. Msi2lmp will be modified again to make use of this feature. NE-1 will read force field files to maintain up-to-date atom types.
Conclusions
Figure 11. Six point glassy Tg least squares.
study could be done to understand the physics in this range and determine a better heating rate and equilibration time. The simulated bulk modulus overshoots by less than 3%. Using the average Tg value calculated (433 K) from the three least squares Tg values determined from Figures 10, 11, and 12, the average Tg value undershoots the experimental value (436 K) by less than 1%. The experimental value was determined by Sindt to be 436 K.[7] The initial density, bulk modulus, and glass transition temperatures are near the experimental values providing a high-comfort level for the MDS software in its alpha-level release state.
On Going Software Development Work on the MDS software is ongoing. The software is under configuration at source forge, http://sourceforge.net/projects/
OSS can be effectively integrated through software modifications, thus, providing new capabilities that the separate applications did not have in their original versions. The MDS software developed in this research, in conjunction with LAMMPS, successfully provided the necessary base geometry for calculating the DGEBA/IPD system densities as a function of temperature, the bulk modulus, and the glass transition temperature to within a few percent of experimentally determined values. MDS, thanks to the robustness of packmol, generated a complex initial cell containing intertwined oligomers. Finally, MDS allowed the user to generate the LGIF for the initial cell in hours, instead of months.
Acknowledgments Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors, and do not necessarily reflect the views of the National Science Foundation. We wish to thank The University of Denver for providing an excellent research atmosphere and for additional funds for Mr. Allen’s graduate teaching assistant appointment. Mr. Allen also wishes to thank The Boeing Company for its support of his masters’ work in materials science, making this article possible, and Dr. Leora Peltz for her early guidance. Keywords: molecular dynamics large-scale atomic/molecular massively parallel simulator software integration open source
How to cite this article: B. M. Allen, K. Paul, M. Predecki, M. Kumosa. J. Comput. Chem. 2014, 35, 756–764. DOI: 10.1002/ jcc.23537
Figure 12. Four point glassy Tg least squares.
[1] S. Plimpton, J. Comput. Phys. 1995, 117, 1. [2] S. J. Plimpton, P. Crozier, A. Thompson, LAMMPS Users Manual, 2013, Available at: http://lammps.sandia.gov/doc/Manual.html, Accessed on December 19, 2009. [3] NIH Center for Macromolecular Modeling and Bioinformatics, VMD Visual Molecular Dynamics, Available at: http://www-s.ks.uiuc.edu/ Research/vmd, Accessed on March 30, 2013.
Journal of Computational Chemistry 2014, 35, 756–764
763
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
[4] C. Wu, W. Xu, Polymer 2006, 47, 6004. [5] M. Simms, Nanorex, Inc., NanoEngineer-1, Available at: http://www. nanoengineer-1.net, Accessed on June 12, 2012. [6] L. Martinez, L. R. Andrade, E. G. Birgin, J. M. Martinez, J. Comput. Chem. 2009, 30, 2157. [7] O. Sindt, J. Perez, J. F. Gerard, Polymer 1996, 37, 2989. [8] J. L. Tack, D. M. Ford, J. Mol. Graphics Model. 2008, 26, 1269. [9] Accelrys, Materials Studio, Available at: http://accelrys.com, Accessed on January 2, 2012. [10] J. R. Maple, U. Dinur, A. T. Hagler, Proc. Natl. Acad. Sci. USA 1998, 85, 5350. [11] H. Sun, J. Phys. Chem. B 1998, 102, 7338.
764
Journal of Computational Chemistry 2014, 35, 756–764
[12] Scooter Software, Inc., Beyond Compare, Available at: http://www. scootersoftware.com/moreinfo.php, Accessed on March 15, 2011. [13] J. Choe, B. Kim, Bull. Korean Chem. Soc. 2000, 21, 419. [14] H. Horstermann, R. Hentschke, M. Amkreutz, M. Hoffmann, M. WirtsRutters, J. Phys. Chem. B 2010, 114, 17013. [15] A. Shokuhfar, B. Arab, J. Mol. Model. 2013, 19, 3719.
Received: 1 November 2013 Revised: 27 December 2013 Accepted: 1 January 2014 Published online on 4 February 2014
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
Integrating Open-Source Software Applications to Build Molecular Dynamics Systems Bruce M. Allen,* Paul K. Predecki, and Maciej Kumosa Three open-source applications, NanoEngineer-1, packmol, and mis2lmp are integrated using an open-source file format to quickly create molecular dynamics (MD) cells for simulation. The three software applications collectively make up the open-source software (OSS) suite known as MD Studio (MDS). The software is validated through software engineering practices and is verified through simulation of the diglycidyl ether of bisphenol-a and isophorone diamine (DGEBA/IPD) system. Multiple simulations are run using the MDS software to create
MD cells, and the data generated are used to calculate density, bulk modulus, and glass transition temperature of the DGEBA/ IPD system. Simulation results compare well with published experimental and numerical results. The MDS software prototype confirms that OSS applications can be analyzed against real-world research requirements and integrated to create a C 2014 Wiley Periodicals, Inc. new capability. V
Introduction
calculating physical properties of an epoxy material. This was achieved by calculating the density of the diglycidyl ether of bisphenol-a and isophorone diamine (DGEBA/IPD) epoxy system over a large temperature range allowing the simulated epoxy to transition from the glassy to viscous state. The glass transition temperature and the bulk modulus of the simulated epoxy were also calculated. The DGEBA/IPD system was selected for simulation due to its complexity and the fact that this system had been previously investigated providing experimental and numerical data for comparison.[7] The DGEBA and IPD structures are shown in Figures 1 and 2, respectively. The actual epoxy material consists of a three-dimensional network of polymer and cure molecules intertwined with other networks of itself creating a crosslinked and intertwined amorphous system. For simplicity, oligomers composed of nine DGEBA and four IPD molecules were used in the simulation in an attempt to model the epoxy material. Each oligomer was 75% intracrosslinked within itself as shown in Figure 3. A total of five oligomers were used to create a cell, which provided intertwining to the cell. Tack and Ford modeled and simulated a similar system consisting of DGEBF and DETDA using the Amorphous Cell in Materials Studio by Accelrys.[8,9] Modeling a similar system using MDS was important due to the intertwining of oligomers provided by packmol.[6] The Consistent Family of Force Fields (CFF91) was used for all simulations.[10] Figure 4 illustrates through a flowchart the three components of the MDS software applications, shown in yellow, and the LAMMPS application shown in green. The integration between the three applications was accomplished with
The goal of this work was to integrate existing open-source software (OSS) applications for the purpose of quickly creating an initial molecular dynamics (MD) cell, also known as a system of molecules or a system for simulation using the largescale atomic/molecular massively parallel simulator (LAMMPS).[1] LAMMPS and the visual MD (VMD) software are used to run MD simulations and to view atom trajectories.[2,3] LAMMPS itself is an object toolkit and can be built as a library and called as part of a larger software application. Colleagues reported initial epoxy system creation times of up to 6 months due to the complexity of specifying single system geometry and associated force field data for LAMMPS. Also, ring spearing and manually fixing systems generated with commercial software made finding alternative system creation approaches attractive.[4] A search for chemistry or nanotechnology computer-aided design (CAD) software yielded NanoEngineer-1 (NE-1), while a review of tools bundled with LAMMPS provided msi2lmp.[5] The initial project direction was to integrate the two software applications, but NE-1 lacked a MD cell or cell concept and drawing medium to large geometries (at least several thousand atoms) would take significant time. Later, while attending the LAMMPS Workshop in August 2011, packmol was brought to our attention, which had the necessary cell concept explicitly designed into it.[6] Packmol was determined to be the missing software application needed to complete the OSS development. The three integrated applications, NE-1, packmol, and msi2lmp are now known collectively as MD Studio (MDS). Initial development requirements were based on the need for application source code and a development environment to support software modifications and test. A lengthy comparison of many OSS applications was not attempted.
Simulated Epoxy System The MDS software integration was tested in this research by following standard software engineering practices and by 756
Journal of Computational Chemistry 2014, 35, 756–764
DOI: 10.1002/jcc.23537
B. M. Allen, P. K. Predecki, M. Kumosa University of Denver, Daniel Felix Richie School of Engineering & Computer Science, Mechanical & Materials Engineering, 2390 South York Street, Denver, Colorado 80210 E-mail: [email protected] Contract grant sponsor: National Science Foundation [(NSF) Grant Opportunities for Academic Liaison with Industry (GOALI)]; contract grant number: #CMMI-1232520 C 2014 Wiley Periodicals, Inc. V
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
Figure 1. Diglycidyl ether of bisphenol a (DGEBA) epoxy monomer resin.
Figure 2. Isophorone diamine (IPD) crosslinker (curing agent).
the development of the enhanced molecular machine part (EMMP) file. The EMMP file was developed for sharing data among the three applications, and is based on the existing molecular machine part (MMP) file. The LAMMPS geometry input file (LGIF) is the final output from the MDS software. The LGIF and the LAMMPS Command Input File were provided to the LAMMPS software, and the LAMMPS output trajectory file was generated and used to calculate the final simulation results (e.g., Tg, density, volume, temperature) for the epoxy system of interest. The trajectory files generated by LAMMPS were viewed using the VMD application.
Software Use Following the flowchart in Figure 4, it took approximately 2 h to draw a single oligomer. Single DGEBA and IPD molecules were drawn and then copied several times. Subsequently, the molecules were attached to one another to form the oligomer.
Figure 3. Oligomer structure.
SOFTWARE NEWS AND UPDATES
Care was taken to not introduce high energies into the starting molecules or when attaching them to form the oligomer. The NE-1 molecular mechanics energy minimzer (MMEM) was executed periodically to lower the molecules’ and oligomers’ potential energies. The existing NE-1 “combine” feature was used to group the molecules in the oligomer into a single molecule. Grouping the attached molecules into a single molecule (oligomer) was necessary because the EMMP file parser added to packmol, expected a single molecule per file. The final EMMP file created was a template for packmol containing the oligomer. The initial cell of five oligomers is shown in Figure 5. Running packmol is an iterative process; therefore, an initial guess of the initial cell size and shape was made based on the size of the oligomer and the number of oligomers. Packmol was run several times, and the volume of the cell was reduced until no solution was obtained. This previous successful execution provided the minimum volume for the initial cell. Packmol took no more than 1 h to generate the initial cell. Finally, the LGIF was created by msi2lmp in a few seconds. The DGEBA/ IPD system required no more than 3 or 4 h to generate the initial cell. The CFF91 force field file required maintenance to add generic values for missing parameters reported by msi2lmp. Maintenance times are outside the execution of the MDS software, but do have an obvious affect on the simulation results.
Software Features and Modifications NE-1 is CAD software created by Nanorex for molecular-level CAD.[5] The source code consists of Python and C. The WingIDE integrated development environment and debugger was used to test Python software changes. NE-1 is capable of modeling and simulating nanomachines, DNA, carbon nanotubes, and many other systems both organic and inorganic. NE-1 stores complex geometry in its MMP file format. The MDS software uses the MMP file as a molecule template. It provides the positions of atoms and bonding topology from which bond lengths, angles, and dihedral angles can be calculated. The MMP file format was expanded in this research to include atom types, thus creating the EMMP file, and was used for building the DGEBA/IPD initial cell. The original NE-1 uses a force field, NanoDynamics-1 (ND-1), which is element based and cannot be used for MD simulations in LAMMPS. ND-1 provides a single atom type for all chemical environments, whereas most force fields used with LAMMPS provide several per atom. Therefore, the modifications to the MMP file were critical. Journal of Computational Chemistry 2014, 35, 756–764
757
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
Figure 4. Flowchart for MDS software.
Modifications to the NE-1 graphical user interface (GUI) were made to support user access and selection of atom types while drawing molecules. Figure 6 shows the force field chooser on the left hand side of the NE-1 GUI. The atom chooser shows oxygen highlighted in blue, and the generic
oxygen atom type is highlighted in the force field chooser. Notice the CFF91 force field is selected in Figure 6. On the drawing canvas, an oxygen atom is highlighted in yellow, to allow the text and numeric atom types to be displayed. In the interest of completing the prototype quickly, atom types are given numerical representations to preclude token parsing issues and allow every element to have a maximum of 100 atom types. The ToolTips dialog is displayed in Figure 7. The dialog allows the user to select whether numerical or character atom types are displayed when atoms are highlighted on the canvas as presented in Figure 6. The user can choose another force field from the drop down box in the force field chooser. Clicking on the CFF91 title button displays another force field to be selected. The COMPASS atom types, not the associated proprietary force field parameters, were integrated into the NE1 force field chooser as a proof of concept, demonstrating how multiple force fields can be integrated into the application.[11] Packmol was developed by Dr. Leandro Martinez at IMECCUNICAMP.[6] Packmol was written in FORTRAN and the GNU FORTRAN compiler and the dynamic data display debugger were used in the CYGWIN environment under Windows 7 to make and test source code changes. The packmol application reads one or more molecule templates and fills a user defined volume with a user specified number of each molecule. Packmol allows the user to generate many unique cells by specifying a unique pseudorandom number seed. This allows the user to create and later recreate a cell based on initial system configuration information.[6] The original packmol software could not read or write MMP files, and was, therefore, modified to read and write the EMMP file. An example of how packmol functions is as follows: a cube containing 1,000 DGEBA molecules and 500 IPD molecules is created using a single-molecule template for DGEBA and single-molecule template for IPD. Given the size of the DGEBA and IPD molecules and the packing constraints, intertwining of molecules is not expected. This is not the case for the oligomers used in this research.
Figure 5. Initial cell: Five oligomers. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
758
Journal of Computational Chemistry 2014, 35, 756–764
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
SOFTWARE NEWS AND UPDATES
Figure 6. New software feature allows choosing the force field. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Figure 7. New software feature allows display of numeric atom types. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
The oligomers used in this work were unexpectedly intertwined, which was an added benefit to cell creation. For the cell shown in Figure 5, packmol calculated a cubedshaped initial cell with ˚. a side length of 55 A The initial cell was then reduced to the equilibrated final cubedshaped geometry, 30.5 ˚ per side at 1 atm and A 298 K using LAMMPS. The final cell in Figure 8 is rendered in unwrapped atom coordinates.[2] The simulation was run using periodic boundary conditions (PBC) which wrap atoms that leave the boundaries of the cell into positions directly on the opposite side of the cell. VMD was used to display the final cell. VMD attempts to display elements using unique colors. Figure 8 illustrates the point that the final cell is considerably denser than the initial cell seen in Figure 5. Packmol optimized the packing of the molecules by calculating a minimum cost function value. The minimum cost is associated with molecules being very close to one another and providing a minimum system volume, while minimizing the repulsion forces between atoms of different molecules. The cost function contained a user defined constraint for maintaining
Journal of Computational Chemistry 2014, 35, 756–764
759
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
Table 1. MDS requirements map. Requirement Experimental density Crosslinked and intertwined cell geometry/moduli Minimize repulsive forces Force field support Calculate atom partial charges Draw molecules Manually type atoms Modify and add atom types Generate LGIF Uniquely identify molecules
NE-1
packmol
Msi2lmp
LAMMPS
X X
X X
X
X
X
X X
X X
X X
X X X X X
The density requirement determines the free volume in the epoxy cell. Crosslinking and intertwining of polymer chains determine moduli of the material and density. The remaining requirements were cited above as software modifications to the three software applications. Table 2 lists the specific source code changes by application and programming language. Figure 8. Final cell: five oligomers after simulated annealing. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
a minimum distance between atoms of different molecules ˚ . When atoms from different molewith an initial value of 2 A cules were too close to one another, high repulsive forces occurred, resulting in atoms moving a large distance in a single simulation time step. LAMMPS terminates under such conditions due to bad dynamics. LAMMPS supports identification of individual molecules to calculate mean-squared displacement (MSD) during simulation. Packmol was modified again to uniquely identify molecules in the EMMP file allowing MSD calculations for small molecules in the cell. The molecule ID number generated in packmol was passed to msi2lmp for inclusion in the LGIF. The msi2lmp software was created by Dr. John Carpenter at Cray Research.[2] Msi2lmp generated a LGIF from the Cartesian coordinate file/molecular data file (CAR/MDF) file pair, which is one file format supported by Materials Studio, a commercial software application.[9] Msi2lmp was written in the C language, and the Microsoft Visual Studio tool was used to modify and test its source code modifications. The msi2lmp software enumerates all unique bond and nonbonded force relationships among atoms of a molecule and between many molecules composing a system. Msi2lmp also populates the force field parameters from a user supplied force field file into the LGIF. The epoxy material simulation required atom partial charges for Coulombic force calculations in LAMMPS. The original msi2lmp transferred the atom partial charges from the (CAR/MDF) file pair directly to the LGIF when processing files created by Materials Studio. The modified msi2lmp software, in this research, calculates the atom partial charges using the method defined by Sun after reading atom data from the EMMP file.[11] The LAMMPS “full atom” model is written to the LGIF, which is required to support the atom partial charges and the molecule ID created by packmol.[2] The above software modifications reflect the new requirements needed to integrate the three applications. Table 1 summarizes the MDS requirements. 760
Journal of Computational Chemistry 2014, 35, 756–764
Software Verification The MDS software developed in this study was subsequently verified by using standard software test methods of white and black box testing, requiring the use of debuggers to examine data values and pathways taken in the modified software. The MDS software was parallel tested using a single DGEBA/IPD oligomer to generate two LGIF files using two separate paths through the modified msi2lmp software. Parallel software testing was performed by converting the EMMP file into a protein data bank (PDB) file, which was then transformed into a CAR/MDF file pair. The CAR/MDF file part was subsequently converted to a LGIF using the existing pathway in the msi2lmp application for string-based atom types. The EMMP file was rewritten into a LGIF using a different pathway in the msi2lmp application that understood numeric atom types. The Beyond Compare software tool was used to compare the two LGIF files for differences.[12] Figure 9 shows the process for the above parallel testing. The two LGIF files were nearly identical, thus showing that the new MDS software created the correct file for simulation. Subtle differences in the two LGIF files were traced back to how NE-1 created errors in the PDB file. The MD verification simulations were executed using the major simulation parameters listed in Table 3.
Table 2. Software modifications. Software modifications
NE-1
NE-1 EMMP file parser (Python) Packmol EMMP file parser (FORTRAN) Msi2lmp EMMP file parser(C) GUI added force field chooser (Python) GUI changes to tooltips dialog (Python) Passes atom types Calculate atom partial charges (C)
packmol
Msi2lmp
X X X X X X
WWW.CHEMISTRYVIEWS.COM
X X
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
Table 4. Single round of simulated annealing at 1 atm. Step 1 2 3 4
Start temperature (K)
End temperature(K)
Duration(ps)
298 298 600 600
298 600 600 298
50 50 200 50
Begin temperature (K)
End temperature(K)
Duration (ps)
298.0 298.0 328.2 328.2 358.4 358.4 388.8 388.8 419.0 419.0 449.2 449.2 479.4 479.4 509.6 509.6 539.6 539.6 570.0 570.0 600.2
298.0 328.2 328.2 358.4 358.4 388.8 388.8 419.0 419.0 449.2 449.2 479.4 479.4 509.6 509.6 539.6 539.6 570.0 570.0 600.2 600.2
10 100 10 100 10 100 10 100 10 100 10 100 10 100 10 100 10 100 10 100 10
Table 5. LAMMPS NPT ensemble. Step
Figure 9. Parallel testing msi2lmp. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
PBC were used to maintain a constant number of atoms over the simulation duration and to isolate the computations from surface effects. The Van der Waals and Coulombic cutoffs ˚ to support the LAMMPS full atom were set to 9.5 and 10.0 A model. The commonly used in organic system simulations, 1 fs time-step, was assumed.[13]
Simulation Methods and Results The following simulation steps were performed to generate the final cell: MMEM, quench (velocity rescaling), and simulated annealing (SA). Table 4 presents the SA details for each of the five rounds. System equilibrium, after SA, was accomplished by running a series of NPT ensembles, where NPT represents holding the number of moles (N), pressure (P) and temperature for the system constant while allowing the system volume to change over time. The ensembles are listed in Table 5. The bulk modulus was calculated by linearly increasing the system pressure from 1 to 5001 atm at 298 K using the NPT ensemble in Table 6. Table 7 lists the calculated density values based on cell dimensions and temperature. The cell dimensions were
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
collected every 0.1 ps, and the calculated density at that temperature was the average of 1,000 data points. The simulated density value was verified by comparing it to the experimentally determined value from the literature for the same resin system tested at ambient conditions (i.e., 298 K and 1 atm).[10] The agreement was within 7% as shown in Table 8. This clearly indicates that MDS provided a high-quality estimate of the resin’s density. The cell dimensions and the bulk modulus of the resin at ambient temperature are presented in Tables 9 and 10, respectively. In addition, a comparison is made in Table 10 between the experimental and simulated values of the modulus. Similar to the density predictions, the error in the modulus prediction is small and was found to be less than 3%. MDS was further verified by determining the glass transition temperature of the resin of interest. Figure 10 presents a graph of the average density data plotted, in black diamonds, against the average temperatures from Table 7. The line in black is the least squares fit of the first five temperatures vs. density values.
Table 3. LAMMPS simulation parameters. Parameter PBC K-space Solver Van der Waals cuffoff Coulombic force cutoff Time step
Value Default PPPM 9.5 A˚ 10.0 A˚ 1 fs
Table 6. Bulk modulus data collection at 298 K NPT ensemble. Step 1 2 3
Begin pressure (atm)
End pressure (atm)
Duration (ps)
1.0 1.0 5001.0
1.0 5001.0 5001.0
100 100 100
Journal of Computational Chemistry 2014, 35, 756–764
761
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
Table 7. Simulated densities by NPT ensemble. Temperature (K)
Lx (A˚)
˚) Ly (A
˚) Lz (A
Volume (cc) 3 1.0e20
Density (gm/cc)
297.55 328.02 358.12 388.63 418.66 449.02 479.55 509.85 539.95 570.57 599.86
30.75 30.74 30.78 30.83 30.79 30.93 31.07 31.39 31.32 31.46 31.56
30.76 30.76 30.79 30.84 30.80 30.94 31.08 31.14 31.34 31.47 31.57
30.77 30.77 30.80 30.85 30.81 30.95 31.09 30.42 31.35 31.48 31.58
2.8954 2.9106 2.9197 2.9344 2.9222 2.9634 3.0042 3.023 3.0783 3.1164 3.14645
1.056 1.0510 1.0458 1.0427 1.0470 1.0325 1.0185 1.0122 0.99404 0.98187 0.9725
The remaining points were used to generate the least squares fit of the remaining six temperatures vs. density values as indicated by the gray dashed line. The first five points were selected to represent the glassy material prior to Tg, whereas the remaining points represent the viscous material. Tg was also estimated by selecting the first six and four glassy points as shown in Figures 11 and 12. Taking the average value of three determined Tgs using the first four, five, and six glassy data points, the average Tg was predicted to be 433 K.
Discussion Particle–particle–particle-mesh (PPPM) performs nonbonded calculations faster than using the Ewald method in a multicore environment. It fails sooner in a multicore environment, due to missing atoms compared to the Ewald method, making PPPM preferable.[2] The results shown in Table 7 took approximately 3 days on a 2.93 GHz 4-core i7 Intel processor using the PPPM option, whereas the same simulation required approximately 5 days using the Ewald option. Using a combination of MMEM, quench, and SA allowed cell conformations to reach a density near the experimental value at ambient conditions. It took approximately 36 h to achieve the ambient density shown in Table 7. Simulation of a system of small molecules would not have required such a long equilibration time. The oligomer size and perhaps the intertwining are responsible for the long preparation time. The oligomers were used to achieve intracrosslinking in the epoxy, due to previous attempts to manually intercrosslink oligomers led to longer equilibration times. The longer times were due to the introduc-
tion of nonequilibrium geometry, when constructing the crosslinks using NE-1. The calculated system density undershoots by 7% compared to the experimental value.[10] Polymer density is also a function of the degree of crosslinking.[14] Shokuhfar and Arab calculated DGEBA/DETA densities for crosslinking densities from 5 to 81%. The experimental values for 16, 25, and 37% are 1.16, 1.19, and 1.13 gm/cc, respectively, whereas the simulated densities are 1.1, 1.11, and 1.12 gm/cc, respectively.[15] In general, simulations undershoot the experimental values. Epoxy materials absorb water vapor, increasing density. There appears to be enough uncertainty in the experimental and simulation methods making determination of monotonically increasing densities with increasing degree of crosslinking difficult. The fifth density in Figure 10 is higher than expected and was most likely due to the heating rate and equilibration time used in elevating the system temperature from 388 to 418 K. The system is transitioning from glassy and viscous, and further
Table 10. Simulated bulk modulus.
D (MPa)
DV (cc)
V0 (cc)
Bulk modulus (GPa)
Experiential bulk modulus (GPa)
506.6
2.90e-21
2.94e-20
5.13
5.0110
Table 8. Simulated vs. experimental density at 1 atm and 298 K. Simulated density (gm/cc)
Experimental density (gm/cc)
Percent difference
1.056 6 0.02
1.13110
761
Table 9. Simulated bulk modulus raw data.
762
Pressure (atm)
˚) Lx (A
Ly (A˚)
˚) Lz (A
Volume (cc)
1.0 5001.0
30.86 24.76
30.87 29.82
30.88 35.92
2.94e-20 2.65e-20
Journal of Computational Chemistry 2014, 35, 756–764
Figure 10. Five point glassy Tg least squares.
WWW.CHEMISTRYVIEWS.COM
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
moleculardynami/, and is chronicled on our blog at http:// moleculardynamicsstudio.blogspot.com/. Several software enhancements are in progress including static and dynamic crosslinking to build better cells, automatic atom typing, support for all open-source atomistic and united atom force fields, and integration of packmol and msi2lmp applications into NE-1 for ease of use. Better crosslinking will result in ease of use and speed up preparing the final cell for simulation. Automatic atom typing will free the user to manually type exceptional cases, where the software is unable to determine the correct atom type. The CFF force fields usually include an “auto equivalence” table to provide missing force field parameters. Msi2lmp will be modified again to make use of this feature. NE-1 will read force field files to maintain up-to-date atom types.
Conclusions
Figure 11. Six point glassy Tg least squares.
study could be done to understand the physics in this range and determine a better heating rate and equilibration time. The simulated bulk modulus overshoots by less than 3%. Using the average Tg value calculated (433 K) from the three least squares Tg values determined from Figures 10, 11, and 12, the average Tg value undershoots the experimental value (436 K) by less than 1%. The experimental value was determined by Sindt to be 436 K.[7] The initial density, bulk modulus, and glass transition temperatures are near the experimental values providing a high-comfort level for the MDS software in its alpha-level release state.
On Going Software Development Work on the MDS software is ongoing. The software is under configuration at source forge, http://sourceforge.net/projects/
OSS can be effectively integrated through software modifications, thus, providing new capabilities that the separate applications did not have in their original versions. The MDS software developed in this research, in conjunction with LAMMPS, successfully provided the necessary base geometry for calculating the DGEBA/IPD system densities as a function of temperature, the bulk modulus, and the glass transition temperature to within a few percent of experimentally determined values. MDS, thanks to the robustness of packmol, generated a complex initial cell containing intertwined oligomers. Finally, MDS allowed the user to generate the LGIF for the initial cell in hours, instead of months.
Acknowledgments Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors, and do not necessarily reflect the views of the National Science Foundation. We wish to thank The University of Denver for providing an excellent research atmosphere and for additional funds for Mr. Allen’s graduate teaching assistant appointment. Mr. Allen also wishes to thank The Boeing Company for its support of his masters’ work in materials science, making this article possible, and Dr. Leora Peltz for her early guidance. Keywords: molecular dynamics large-scale atomic/molecular massively parallel simulator software integration open source
How to cite this article: B. M. Allen, K. Paul, M. Predecki, M. Kumosa. J. Comput. Chem. 2014, 35, 756–764. DOI: 10.1002/ jcc.23537
Figure 12. Four point glassy Tg least squares.
[1] S. Plimpton, J. Comput. Phys. 1995, 117, 1. [2] S. J. Plimpton, P. Crozier, A. Thompson, LAMMPS Users Manual, 2013, Available at: http://lammps.sandia.gov/doc/Manual.html, Accessed on December 19, 2009. [3] NIH Center for Macromolecular Modeling and Bioinformatics, VMD Visual Molecular Dynamics, Available at: http://www-s.ks.uiuc.edu/ Research/vmd, Accessed on March 30, 2013.
Journal of Computational Chemistry 2014, 35, 756–764
763
SOFTWARE NEWS AND UPDATES
WWW.C-CHEM.ORG
[4] C. Wu, W. Xu, Polymer 2006, 47, 6004. [5] M. Simms, Nanorex, Inc., NanoEngineer-1, Available at: http://www. nanoengineer-1.net, Accessed on June 12, 2012. [6] L. Martinez, L. R. Andrade, E. G. Birgin, J. M. Martinez, J. Comput. Chem. 2009, 30, 2157. [7] O. Sindt, J. Perez, J. F. Gerard, Polymer 1996, 37, 2989. [8] J. L. Tack, D. M. Ford, J. Mol. Graphics Model. 2008, 26, 1269. [9] Accelrys, Materials Studio, Available at: http://accelrys.com, Accessed on January 2, 2012. [10] J. R. Maple, U. Dinur, A. T. Hagler, Proc. Natl. Acad. Sci. USA 1998, 85, 5350. [11] H. Sun, J. Phys. Chem. B 1998, 102, 7338.
764
Journal of Computational Chemistry 2014, 35, 756–764
[12] Scooter Software, Inc., Beyond Compare, Available at: http://www. scootersoftware.com/moreinfo.php, Accessed on March 15, 2011. [13] J. Choe, B. Kim, Bull. Korean Chem. Soc. 2000, 21, 419. [14] H. Horstermann, R. Hentschke, M. Amkreutz, M. Hoffmann, M. WirtsRutters, J. Phys. Chem. B 2010, 114, 17013. [15] A. Shokuhfar, B. Arab, J. Mol. Model. 2013, 19, 3719.
Received: 1 November 2013 Revised: 27 December 2013 Accepted: 1 January 2014 Published online on 4 February 2014
WWW.CHEMISTRYVIEWS.COM
