Cross-linked Polymers: Part II - Dynamic Method I

Welcome to part two of cross-linked polymers. I will cover a dynamic cross-linking approach where LAMMPS will put resin and curing agent molecules into motion and molecules will be bonded together external to LAMMPS. This approach is based on LAMMPS as a callable library from a larger cross-linker application. If anyone has invented this wheel why not at least give me some feedback. While reading the literature on molecular dynamics, it is easy to find descriptions of cross-linking. Many authors speak to creating some automation code and seem to never publish it. I suppose in order to be a molecular dynamic scientist one must fashion his or her own cross-linker  - just as a Jedi must fashion his or her own light saber.  Another piece of vaporware I am interested in building will do the following to create an amorphous polymer matrix:

• The cross-linker application will defined an internal geometry data structure for all molecules taking part in the simulation. The data structure will consist of sub-structures defining molecules, atom types, individual atoms, bonds, angles and dihedral angles. All of this information will change over time as bonds are broken and others are created. As the cross-linking occurs over time the number of molecules will decrease! 
• The "start points" previously defined in part one, to attach molecules, will be used in this new software for the same purpose. A DGEBA molecule will have two start points and an IPD molecule will have four start points. The user will define the probability of when an IPD molecule is attached to a DGEBA molecule and the same information for DGEBA molecule for an IPD or another DGEBA molecule. Keep in mind that your molecules my be different - just create a template for each one. 
• A radius of reaction will be defined by the user for every start point. Only start points within the radius of reaction can be considered for cross-linking. Perhaps you flip a fair or unfair coin to determine if the cross-linking occurs or not. The software will be flexible enough to support this.
• A number of LAMMPS simulation steps will be defined by the user so that LAMMPS can be executed for this number of steps. The number of steps and the time duration could be based on your knowledge of the system kinetics.  
• A maximum number of LAMMPS executions will be defined so that the cross-linker software does not execute forever. 
• A desired cross-linking percent will be defined by the user to tell the cross-linker how far to execute.
• A user defined convergence value will be used to test cross-linking percentage so the cross-linker can terminate in the event the desired cross-linking could not be obtained.
• If specified, the cross-linker will tell LAMMPS to run a molecular mechanics energy minimizer, MMEM, to help take some of the stress out of the geometry created by the new bonds. Perhaps, prior to running the MMEM the software could be instructed to run the structure cleaner. I would like to see quantify how useful the structure cleaner could be in reducing times to achieve low potential energy values.
• The cross-linker will call LAMMPS to stir up the soup of resin and curing agent molecules and get the system away from bad dynamics.
• The cross-linker will run LAMMPS for some user defined number of simulation steps and return to the cross-linker so it can test molecule start points against one another for possible cross-linking.
• A new bond will be formed if the radius of reaction criteria is met and the number of molecules in the cross-linker will decrease retyping atoms in the area of the new bond(s) will occur. 
• Some atoms will be deleted due to the cross-linking and will be removed from the simulation.
• Termination criteria will be tested.
• LAMMPS will be executed again from inside the cross-linker program.

The above approach to "chemical reaction" will result in an amorphous cross-linked polymer. I intend to compare it to an approach using the ReaxFF. Additional infrastructure software will have to be modified to accomplish the goal of amorphous cross-linked polymers. One interesting feature of PACKMOL is that it packs the molecules given to it in a serial fashion. The first molecule is packed then the second and so on. This leads to a segregated MD cell. PACKMOL requires an enhancement to distribute all molecules across the entire MD cell volume. This will provide the best mixed environment for cross-linking. One work around is to create a PACKMOL input file that specifies DGEBA and IPD in alternating references for the total number of desired molecule in the MD cell. It you need one hundred of each this works well. 

Another use of PACKMOL would take advantage of creating a fixed molecule structure, such as the cross-linked polymer created above, and then instruct PACMOL to distribute additional molecules into the structure. Fixed molecules are used only to take up space in the MD cell and the new molecules added will be no closer than the tolerance limit provided by the user. This approach will require converting the final cross-linked polymer MD cell into a molecular machine part (MMP) file so PACKMOL can read it. A small molecule could then be packed into a small portion of the MD cell that contains the cross-linked polymer. Water could be placed into a well defined MD cell at desired temperature and density. 

The dynamic cross-linking approach will be prototyped by building software in C++ to do specific tasks. The "knitter" application will read a LAMMPS output file and build the MD cell geometry in MMP file format in order to check for "start points" that meet the reaction radius requirements and "knit" molecules together. Next, a new MMP file containing new geometry, atom types, molecules, etc. will be created. MSI2LMP will be used to create the new LAMMPS Geometry Input File and LAMMPS will be executed again base on the desired number of steps. All the applications involved will be called from a script.

Cross-linked Polymers Part III will examine an idea that could allow PACKMOL to generate cross-linked systems as it packs molecules. PACKMOL uses constraints to put molecules as close together as possible without violating the user defined tolerance - usually 2 angstroms, which stops atoms from different molecules from being too close. In a nutshell, what if a dynamic constraint could be added to PACKMOL so that after adding DGEBA an IPD molecule is added but the "start point" on the DGEBA defines the position for a start point on the "to be added" IPD molecule. 



    

Cleaning House & May the Force be with You!

Just a quick update on NanoEngineer-1. I have installed Windows 7 SP1 just to test NanoEngineer-1 in XP compatibility mode. Sadly, this was a failure. As far as Windows 7 is concerned I still have ditch the service pack. If you have any additional information on this problem please let me know. I also tried modes for Windows 7 and 2000 with no luck.


I briefly mentioned potential energy of the molecules in my last post. NanoEngineer-1 has GROMACS molecular dynamics simulator integrated into it. This would be a great idea for LAMMPS too. You can draw a molecule and start up the simulator and then watch a movie of the molecule under dynamics. The cool part of this is bonds can be stretched to hugh lengths and the simulator does not break the bond during simulation. Instead the molecule experiences a hugh influx of energy and bounces around. How does one stop such behavior? One solution is to minimize the energy of the molecule(s) using the molecular mechanics energy minimizer provided by NE-1! The MMEM uses some algorithm such as steepest decent and calculates a local energy minimum for the system or for specific highlighted molecules. MMEM uses the NanoDynamics-1 force field. If memory serves me right it minimizes on bonds, angles and does not pay attention to dihedral angles. It appears to pay attention to molecules that are not bonded together and perhaps has something like Van der Waals forces. The down side of the MMEM is it takes a long time for large molecules and since you are only guaranteed a local minimum in energy you may be starting from a weak point on the potential energy surface. So waiting a long time and probably not getting what you want is no fun. It's time to go shopping for a better solution.

There are many chemistry applications on planet Earth and one is ChemSketch by ACD/Labs. I used it to draw molecules for school. I noticed that it had a "cleaner" feature. Not having the source code I can't look at how it does this cleaning of the structure so quickly. Knowing something about chemistry I can guess how the cleaner works. Let's face it atoms generally hate each other. A methane molecule has hydrogens spread out such that all five atoms are as far away from one another as possible. Molecules appear to want to maximize the space they occupy. So one idea is to allow atoms to assume ideal positions relative to one another based on hybrid orbitals. Sp1, Sp2 and Sp3 provide simple geometry to be used in placing atoms relative to one another. Of course not all atoms or groups of atoms are treated the same by a force field. Looking at COMPASS family of force fields one sees many equilibrium bond lengths between atoms based on their atom types. So the hybrid orbitals could be a good first step, but using the equilibrium bond lengths, bond angles, and dihedral angles from the desired force field would be even better! Besides a computer has nothing better to do than to recompute (x,y,z) positions for an atom relative to another atom based on the equilibrium values in the force field. Some of this data will be missing but the fall back would be to substitute generic values. 

The molecule template created by NanoEngineer-1 contains cartesian coordinates. There is nothing wrong with this data format, but the force field data is in a format that chemists recognize immediately! The Z-matrix is used to represent a molecule or a system of atoms/molecules by their natural coordinates - bond length, bond angle, dihedral angle and out-of-plane angle. If you have used The Gaussian software you know about the Z-matrix. So in the interest of creating bad puns: in order to clean house you need the force to be with you! 

My next piece of vaporware will be a cleaner software and hopefully be implemented in Python inside NanoEngineer-1. It will have a cleaner button on the GUI and will do the following behind the scenes: first a Z-matrix will be constructed from the cartesian coordinate data and some dummy atoms will  be added to support the Z-matrix. The Z-matrix data structure is devoid of cartesian coordinates. The bond lengths, bond angles, dihedral angles and out-of-plane angles could be calculated and stuffed into the Z-matrix, but why bother. These values will be replaced with the correct data from the desired force field. Every atom carries an atom type from the force field and this allows the look up of the correct values to fill in the Z-matrix. Once completed, the Z-matrix will be transformed into a molecular machine part (MMP) internal representation within NanoEngineer-1. A single pass through the atoms should create the Z-matrix and a second pass should create the MMP internal representation. The effort is linear 2N!

A final pass may be necessary to ensure that all atoms are as far away from their neighbors as possible. Since we are not calculating energy and doing a steepest decent minimization we may need to develop an objective function based on maximizing distance between atoms while maintaining equilibrium bond lengths, etc... The above is the easy part of the problem. The Van der Waals and Coulombic contributions in minimizing the overall system energy appear to be satisfied only by a full MMEM. Perhaps the structure cleaner will serve to provide a better starting point. One obvious use is to clean up manually cross-linked molecules drawn in NanoEngineer-1. I had manually cross-linked a few oligomers to attempt to get higher densities for the DGEBA/IPD system.  The simulations ran much longer due to structure "dirtiness" created by my really long bond lengths. Since the oligomers contain cycles the software must handle such cases and not step through the structure forever! 

Comments and suggestions are invited and encouraged here! 

Cross-linked Polymers: Part I - Static Method

Chemical reactions are fun! Throwing some resin monomers into a beaker and mixing them with a curing agent results in something greater than the sum of their individual molecules. Once it cures you have a material that should have been jammed between two things you wanted glued together or something you should have poured into a mold to make one of those vendetta masks I see on TV these days. All those molecules magically react with one another forming a 3D-network of molecules. monomers form bonds with other monomer molecules and with curing agent molecules.

The world of molecular dynamics (MD) does in some cases support chemical reactions. The LAMMPS MD simulator supports the ReaxFF and sadly I have not devoted much time to it. Ignorance can be a dangerous thing. One could reinvent the wheel and look foolish especially if the wheel produced is not as good as a standard wheel. For what it does ReaxFF is a fantastic tool! It is 10-50 times slower than non-reactive force fields, but is NlogN in performance compared to quantum mechanical (QM) methods which are N^3! This is according work done by Adri van Duin in a power point I read online. Ideally, it would be fantastic to build a cell filled with DGEBA and IPD molecules and specify the ReaxFF force field and let LAMMPS build an initial cell to some desired level of cross-linking. ReaxFF is on my to do list. I am working on less creative and more brute force methods for creating cross-linked systems, which I plan to benchmark against the ReaxFF approach.

My first attempt at cross-linking software is still vaporware. It is based on a static approach. The static approach is based on building molecules one atom at a time using a template to tell the software where to place the next atom relative to the previous atom. Starting points in the cell are randomly generated initially and a molecule to be built per starting point is randomly decided. Atoms are added to each of N molecules in a round-robin fashion giving nearly equal access to each molecule under construction. Each molecule is built in turn and positions on a molecule's structure is identified as a point where another molecule can be attached. If there are many molecules to choose from for this attachment then a random selection is made based on the odds supplied to the software. A DGEBA molecule has two places in its structure where either another DGEBA or an IPD can be attached. If the choice is 50:50 then a coin is flipped to choose one. Each DGEBA molecule contributes two new starting points to a stack of starting points in the software. The situation could become explosive! Every IPD molecule has four start ing points so explosive is an understatement! So what allows the software to terminate? The software uses several methods of termination. The first is to achieve a desired cell density. The second is collisions. Every molecule under construction could  collide with another molecule being built! Which ever molecule that tests for and finds a collision is selected for deletion! This applies to molecules that attempt to leave the cell. The third and fourth tests used to terminate the molecule building process is to track the cell density and upon convergence to any value or after some user defined iteration limit the process terminates.

What is collision detection? In the case of molecules it means that two atoms cannot occupy the same space. It also means that a bond cannot spear a ring structure. An atom should not occupy space between two bonded atoms. The book, "Real-time Collision Detection" by Christer Ericson provided me the knowledge about collision detection. Testing for collisions between points with points, lines (aka bonds)  and planes (aka ring structures) , and lines with planes are being done in computer games and so it can be used in a much slower environment of building cells for simulation. Organic ring structures are represented by one or more unique triangular planes. Again, collision detection causes the current molecule to be erased and the starting point to be replaced with a hydrogen atom. The following  structure shown below shows DGEBA like molecules being built atom-by-atom where the monomers are being assembled into chains simulating a polymerization process. Collision detection is not the Holy Grail. Software is not smart and every new atom, bond and plane needs to be tested for collisions against all the other atoms, bonds and planes!!! This is at least N^2 in effort! Fortunately, there are tricks discovered by smart people that can be employed to this problem. The one I use is to partition my cell into many cubes. Every atom, bond and plane belongs to at least one of these cubes. Testing against the atoms, bonds and planes in the cube where the new atom is placed drastically reduces the number of tests. Also, since the atoms in a molecule template do not collide with other atoms in the template by definition these comparisons can be ignored as well.

"Polymerization" of DGEBA like molecules

Notice in the image above that the monomers are in regular patterns that could possibly be the start of some crystalline or semicrystalline structure. This is a side benefit from the static cross-linking and could possibly result in other test cell geometries beyond amorphous cells. Randomly generating a new trajectory angle for each monomer attached to a chain could result in an amorphous structure assuming the chain does not turn back on itself thus negating a self avoiding walk or SAW. In face the name of the static cross-linking software is SAW.

SAW creates trees of molecules bonded together growing throughout a volume. This process does not create a graph where a graph is a tree with multiple paths to a node and cycles. Each independent tree does not eventually cross-link to another in the cell. This is like a dendritic structure. Of course with enough chains entangled this approach to cross-linking may be good enough? As I pointed out on a previous post I used the oligomer approach to cross-linking for my graduate thesis work. PACKMOL packed the five oligomers into a cell and as an extra bonus entangled the oligomers.

Another important issue is the potential energy of the molecules. If one draws the molecule template in NanoEngineer-1 with longer or shorter than expected bonds or angles that are too wide or narrow it is like storing energy in a spring. NE-1 is used to create the molecule template and the stress added will be placed into every molecule replicated by SAW. This energy is released at the beginning of simulation in LAMMPS and leads to failed simulations because atoms jump across CPUs. Bad dynamics are bad! I will have more to say about energy in a separate blog article soon.

In part two I will cover a dynamic cross-linking approach where LAMMPS will put resin and curing agent molecules into motion and molecules will be bonded together external to LAMMPS. This approach is based on LAMMPS as a callable library from a larger cross-linker. If anyone has invented this wheel why not package up your source code and make a contribution to open source? Comments and suggestions are very welcome on this blog. You can't hurt my feelings because I have worked in the software industry for over 20 years and have a tough skin. Let's work together to build tools to build better students who will go to industry and have some understanding of real-world problems.

Closeup of DGEBA like structures above

Happy New Year 2013!!! PACKMOL, Chunks, Molecules & LAMMPS

I found an embarrassing bug in PACKMOL. Actually, in the modifications I made to PACKMOL. How could this happen? Writing software is not glamorous and often thankless. At the best of times the user should get their work done and at the worst of times the user is ready to fire the programmer! Software is defined as "code, data and documentation maintained in a specific configuration." I have a software engineering book that explains this. Explicit requirements for software are code, data and documentation. The implicit requirement is the former are maintained in a specific configuration. Sadly, I lost configuration! I had fixed this problem a while back and lost the changes. I have been working behind the scenes to build a better infrastructure to prevent this. I will soon host all my source code on Dropbox. I will have a directory structure on Dropbox and point my development tools at Dropbox. Of course this sounds silly to many software engineers. Why didn't I just host my project on github and forget about Source Forge? Why am I not using state of the art development practices? I will change my evil ways once this project attracts additional software developers and collaboration becomes critical. As of this writing, the software is alpha and in need of additional developers who want to own this project.

Back to the bug I found. My modified PACKMOL reads and writes molecular machine part (MMP) files. This is a file format created by the NanoEngineer-1 molecular CAD application. I misplaced memory location of the molecule templates for all molecules after correctly processing the second molecule template. I also created a phantom molecule in the MMP file. I fixed this bug as well. Now the MMP file created by PACKMOL assigns a "chunk" to every molecule. What is a chunk? A chunk is a concept from NanoEngineer-1. It identifies one or more atoms. When you click on a chuck identifier inside the NE-1 CAD program the collection of atoms are highlighted in green. You can now move the entire chunk around on the screen. The funny thing is that the "mol" token is used in the MMP file to define a chunk. Someone must have originally wanted to make a chunk a molecule, but that did not happen completely. NanoEngineer-1 creates chunks willy-nilly. There can be one or more atoms in a chunk. A single molecule drawn in NanoEngineer-1 can contain many chunks as you create atoms as well as copy and paste existing atoms in order to save time drawing a molecule.

When using NanoEngineer-1 the user must enter "Chunks" mode by clicking the chunks button. Upon existing chunks mode one can see the tree of chunks down the left-hand-side of the application. Again, clicking on one or more will highlight the chunks in the canvas showing the chunks in the molecule drawn. Selecting more than one chunk is done by holding down the shift button and clicking on the chunks in the tree. Looking at the icon menubar one sees an icon of three green squares and a black arrow pointing to a single green rectangle. The tooltip tells you this is a chunk combiner tool. If you use this tool to combine chunks into a single chunk then all combined chunks will be highlighted. This is important in using LAMMPS!

LAMMPS has built into it the concept of a group-ID. This is a way to identify a group of atoms. For my purposes a group would be a molecule such as a water molecule floating around in a polymer matrix. LAMMPS also supports a full atom model. The model includes a molecule ID on every atom of a simulation.    The molecule-Id can be used as a group-ID. One can parse LAMMPS output files for specific molecules and use trajectory data to create a path a molecule travels over the course of a simulation. In order to get this molecule data to LAMMPS correctly from NanoEngineer-1 and PACKMOL you need to follow some simple rules:

1. Draw a single molecule per MMP file you create inside NanoEngineer-1. Once you have your single molecule completely drawn then combine the chunks into a single chunk.
2. Let PACKMOL create your initial MD cell from multiple MMP files. The MMP file processing in PACKMOL will uniquely identify all molecules created from the MMP files.
3. MSI2LMP reads and understands the chunks in the MMP file created by PACMOL and will assign unique molecule-IDs to all molecules when it creates the LAMMPS Geometry Input File.

I hope the fixes to PACKMOL and this discussion about chunks and molecules help you in creating successful MD cells. I have seen many downloads to date, but have not received many, three, comments back. Open source is about helping people be successful by making quality software. I look forward to comments, and feature suggestions. PACKMOL should create quality MD cells from many molecules now.