Monday, July 8, 2013
Running LAMMPS to run LAMMPS...
I naively believed one could build a system of molecules outside of LAMMPS and then pop the LAMMPS Geometry Input File into LAMMPS and run an ensemble (e.g. NPT, NVT, etc...) and have answers pop out the other end. I was first confronted with instantaneous pressure values which are alarming if you are new to molecular dynamics (or perhaps ignorant as I was). The pressure values are a topic for another blog entry. Looking at MD cell creation as the context of "running LAMMPS to run LAMMPS" one can say that we are all in the hunt for a free lunch in terms of a minimal total wall-clock time to achieve a correct answer from an MD simulation. The MD cell is the starting point. Build an MD cell close to physical reality and simulation can begin immediately. the devil is in the details.
I have written about several methods to cross-link molecules. The static method (see Cross-linked Polymers: Part I - Static Method) starts with an empty MD cell and attempts to create polymer chains composed of molecules atom-by-atom until the desired density is obtained. The desire is to have many polymer chains intertwined by this process in addition to the cross-linking that will occur when a monomer molecule is completely built and a curing agent molecule is started from its "hotspot" atom. A 3D network of monomers and curing agent molecules is built for each starting point in the MD cell. Since all polymer chains have different attitudes they should intertwine. The missing piece from the static cross-linking puzzle is the final system density - the target objective. Static cross-linking occurs without MD. Therefore, the kinetic energy component is missing. What happens when the system undergoes MD simulation? Hopefully, the short distance between atoms of different molecules does not cause bad dynamics. If the system does not achieve the desired density at state (i.e. temperature and pressure) are there too few degrees of freedom in the system? Is the system non-ergodic? Is it possible to constrain the placement of atoms of different molecules as some safe distance just as packmol does? An atom could be rotated such that it maintains its correct bond length, but this could force a recalculation for the positions of the other two hydrogen atoms in say a methyl group which were placed successfully prior to the atom in question.
Static cross-linking is the Yin of MD cell creation while a dynamic cross-linker (see Cross-linked Polymers: Part II - Dynamic Method I and Dynamic Method II) is the Yang. Dynamic cross-linking starts with an MD cell full of monomers, hopefully distributed to foster cross-linking without the need for significant diffusion of molecules where diffusion requires time. Dynamic cross-linking immediately solves the kinetic energy component by placing molecules at the correct distance apart for nonbonded conditions. Creating a bond between two molecules to partially create a polymer chain will reduce the volume occupied by the polymer since the bonded molecules are closer together. The random attitudes of the monomer and curing agent molecules ensure a polymer chain that follows a random path through the MD cell. Both MD and molecular mechanics energy minimization (MMEM) are required for every group of new cross-links created. Where the static method could recalculate conflicting atom positions using collision detection and rotation of atoms around a common axis the dynamic method avoids this through MMEM. Once cross-links are placed MMEM attempts to equilibrate the system. MMEM is used because there is no analytic function relating geometry and potential field to state (i.e. temperature and pressure). As monomers and curing agent molecules are connected to form polymers, their mobility decreases. The final system density, rate of cross-linking and final percentage of cross-linking is dependent on diffusion of the monomers and curing agent molecules. Could the initial distribution of monomers and curing agent molecules prevent the formation of a system with the desired density in a timely manner? Is the system non-ergodic?
Both approaches create MD cells, but their final geometries will determine the amount of wall-clock time is needed to run LAMMPS for the system to be deemed valid for simulation. Perhaps simulated annealing (SA) is required to achieve a valid MD cell ready for simulation. I currently use SA to build my initial MD cell, which uses packmol with a MD cell full of oligomer rings. The rings are cross-linked monomers and curing agent molecules forming oligomers. Half my simulation time is spent in MD cell preparation via SA. It is used in conjunction with the NPT ensemble to achieve a desired density at state. Both methods discussed above lack an analytical function to relate geometry and potential field to state.
Therefore, running LAMMPS to run LAMMPS is necessary - at least for me since I don't have the magic software to build fully equilibrated systems directly. I'm not complaining. Whether you are laughing or agreeing with my blog entries please drop me a line and let me know your opinion. I am finally starting on the design and code for a dynamic cross-linker today, July 8th, 2013.
Tuesday, May 21, 2013
State of the Union
I suppose one of the main reasons for vaporware is the unintended situations we all find ourselves in the coarse course of daily events. I wanted to update this blog weekly and have had demands on my time I didn't expect - the first being medical. Definitely not life threatening, it continues to be life coloring. It has hampered my progress in my pursuit of knowledge, my research and of course my software development. I had hope to address a cross-linker software by now. That is my objective this summer in addition to running laboratory experiments in support of my molecular dynamics modeling.
I just finished a two quarter class in Electrodynamics! What a rush! I have to admit many Materials Scientists walk around with one eye shut. One eye on mass and the other shut to charge. If anything my small simulations of DGEBA/IPD showed me the power of partial charges on matter and why Van der Waals forces are necessary for good simulations. I now have a much better understanding of polarization of materials in an electric field. I think it cleared up some questions I had about polarization in chemistry too. I suggest to all Materials Science majors to take Electrodynamics. Imagine the extra "dog pile" one could add to a molecular dynamics simulator - not "that" kind of dog pile, but the one where all the dogs jump on you because you have a hotdog in your hand. As the geometry of your molecule changes with time the partial charges could be recalculated. The electric field in the material would then change due to the redistribution of charge. The moving charges give rise to magnetic fields which can induce electric fields and vs versa! Griffiths' Electrodynamics text points out the wildly fluctuating fields inside matter at the microscopic level. How these fields affect diffusion, alignment of molecules, whether the fields help or hurt molecules to clump could affect reactions of water with polymers, and how charged molecules behave inside matter. Also, when materials are in a significant electro-magnetic field perhaps fun things happen to the molecules that are diffusing in the materials.
This summer will provide freedom to work on the following:
I just finished a two quarter class in Electrodynamics! What a rush! I have to admit many Materials Scientists walk around with one eye shut. One eye on mass and the other shut to charge. If anything my small simulations of DGEBA/IPD showed me the power of partial charges on matter and why Van der Waals forces are necessary for good simulations. I now have a much better understanding of polarization of materials in an electric field. I think it cleared up some questions I had about polarization in chemistry too. I suggest to all Materials Science majors to take Electrodynamics. Imagine the extra "dog pile" one could add to a molecular dynamics simulator - not "that" kind of dog pile, but the one where all the dogs jump on you because you have a hotdog in your hand. As the geometry of your molecule changes with time the partial charges could be recalculated. The electric field in the material would then change due to the redistribution of charge. The moving charges give rise to magnetic fields which can induce electric fields and vs versa! Griffiths' Electrodynamics text points out the wildly fluctuating fields inside matter at the microscopic level. How these fields affect diffusion, alignment of molecules, whether the fields help or hurt molecules to clump could affect reactions of water with polymers, and how charged molecules behave inside matter. Also, when materials are in a significant electro-magnetic field perhaps fun things happen to the molecules that are diffusing in the materials.
This summer will provide freedom to work on the following:
- Health - new medication, continue my exercise program
- Research - experiments into a variety of ways to damage siloxane materials and how to model such damage so changes in physical properties can be seen in MD simulations
- Software - build a cross-linker
- Research - how does electrodynamics fit into molecular dynamics simulations
- Research - which force field or potential will work best for simulating siloxane materials such as PDMS
I will present weekly progress on the cross-linker and the technical details. Also, extending existing software into millions of atoms should prove interesting in terms of time efficient algorithms. I strongly believe that MDS could be a great multi-disciplinary teaching project across many fields of study. I look forward to discussing enhancements with all who are interested in enhancing this open source project.
Saturday, February 23, 2013
LAMMPS Wish List
Large-scale Atomic/Molecular Massively Parallel Simulator or LAMMPS is a fantastic research tool! When one looks at the list of features it's amazing that LAMMPS is so fantastic. One feature is the electric field. The following YouTube video shows the power of LAMMPS combined with an electric field capability.
https://www.youtube.com/watch?v=lMBRnxbcRH4
How could LAMMPS need a wish list? Well I do have a wish list for LAMMPS. I'm not being ungrateful, I am wishing and willing to help add the following features to LAMMPS for my research. Features can be added to the proper source tree or to USER features. The following list will be updated as time goes on.
1. Magnetic fields - adding this feature would allow small charged molecules to transport differently in a polymer matrix when in the presences of an electromagnetic field. Water moving in a polymer matrix comes to mind.
2. Ultra-violet radiation damage and cross-linking capability - adding this capability would allow the study of a decomposing polymer matrix, the transport of low molecular weight material created due to damage and cross-linking in polymers. ReaxFF would be a good first starting point. I was thinking of bombarding the MD cell with phantom particles that would interact with atoms to cause bond scission and allowing ReaxFF to form a new bond or reform the broken bond.
I know this is a short entry, but I wanted to get a placeholder for the topics out here. If you have comments I would enjoy hearing them.
https://www.youtube.com/watch?v=lMBRnxbcRH4
How could LAMMPS need a wish list? Well I do have a wish list for LAMMPS. I'm not being ungrateful, I am wishing and willing to help add the following features to LAMMPS for my research. Features can be added to the proper source tree or to USER features. The following list will be updated as time goes on.
1. Magnetic fields - adding this feature would allow small charged molecules to transport differently in a polymer matrix when in the presences of an electromagnetic field. Water moving in a polymer matrix comes to mind.
2. Ultra-violet radiation damage and cross-linking capability - adding this capability would allow the study of a decomposing polymer matrix, the transport of low molecular weight material created due to damage and cross-linking in polymers. ReaxFF would be a good first starting point. I was thinking of bombarding the MD cell with phantom particles that would interact with atoms to cause bond scission and allowing ReaxFF to form a new bond or reform the broken bond.
I know this is a short entry, but I wanted to get a placeholder for the topics out here. If you have comments I would enjoy hearing them.
Cross-linked Polymers: Part III - Dynamic Method II
Welcome to Cross-linked Polymers Part III which will examine an idea that could allow PACKMOL to generate cross-linked systems as it packs molecules. Adding such a capability to PACKMOL would be a great leap forward in building pseudo-random geometries that are in fact cross-linked!!! But first, I want to expand on an idea of fixed atoms from last time. Previously, I stated using a cross-linked MD cell in the fixed mode PACKMOL provides thus locking the polymer matrix in place and then adding water or other molecules into the space remaining in the polymer matrix. After having such a thought I only recently asked why couldn't one distribute monomers into the MD cell using a large tolerance between atoms of different molecules using PACKMOL and then run a second time with the monomers fixed and then have it pack cross-linking molecules into the remaining empty space? Being no free lunch, we must use the first output MMP file generated to create a PACKMOL input file specifying the coordinate of each monomer's baricenter (not the center of mass) and it's Euler angles in order to rotate it into the correct attitude to match the MMP file. This is not to horrible to generate. For thousands of monomers there will be the following entry in the PACKMOL input file:
structure monomer.mmp
number 1
center
fixed x,y,z,a,b,g
end structure
Obviously, either a new piece of vaporware is needed to generate the fixed monomer portions of a second PACKMOL input file or PACKMOL could be extended to make multiple passes in order to build up an MD cell.
In this case a second pass would recreate the monomers in fixed geometry and then the cross-linker molecules would be packed into spaces around them. An important question to ask about dynamic method one is where does the well mixed MD cell consisting of monomers and cross-linking molecules come from? As stated before the current PACKMOL packs the first molecule template followed by the second and so on. This segregates the molecules. If polymerization continues using dynamic method one we will eventually have a system that will not allow the diffusion of monomers or cross-linking molecules. Building a PACKMOL input file of alternating DGEBA and IPD molecules is doable with a Perl or Python script. An now back to our irregularly schedule blog entry...
PACKMOL uses constraints to put molecules as close together as possible without violating the user defined tolerance. 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. Remember, the molecular machine part file (MMP) created in NanoEngineer-1 can contain a "starting point", which I refer to as a "hot spot" If you think of a monomer in a polymer it has two hot spots. They are located at the sites where they react with other molecules. For my system of DGEBA and IPD these are at opposite ends of the monomer. DGEBA has two hot spots - male and female. IPD has four. IPD has two amine groups that have two hydrogens each. One IPD bonds with a maximum of four DGEBA. Therefore, if one wants to create a cross-linked epoxy system at 50 percent cross-linking one would specify only two of the IPD hot spots could be bonded and the number of IPD and DGEBA molecules would need to specified so that excess of either molecule does not occur - unless your system wants to study such a situation.
PACKMOL software changes could be moderate to complex. PACKMOL would pack molecules in a different manner than today. PACKMOL would have to understand the hot spots. When it encounters a hot spot it would use the coordinates of the hot spot to place the other molecule. The molecule to place would be based on the information stored in the hot spot when created in NanoEngineer-1. If there were a 50 percent probability of placing another monomer or a cross-linking molecule at the hot spot PACKMOL would flip a coin. If an atom has multiple hot spots and one has a better chance of reaction then PACKMOL will use this information stored in the hot spot. This new approach to accessing molecules for packing will not affect the packing algorithm. So what does this have to do with building cross-linked MD cells?
PACKMOL currently uses static constraints. One defines the MD cell dimensions and the tolerance between atoms of different molecules. One can select many 3-D MD cell shapes (e.g. spheres, boxes, etc...). Adding hot spots to PACKMOL will add a dynamic constraint. Every monomers would add two dynamic constraints as to where two IPD molecules have to exist. Specifying the placement of the hot spot of molecule B at tolerance distance from hot spot on molecule A sets up the necessary condition for cross-linking the two molecules. Keep in mind that PACKMOL will pack the molecules and maintain the value of the cost function over many attempts and then select the optimal packing configuration as the solution. Keep in mind that only the position of one hot spot is constrained while the attitude of the molecule to be place will vary in the packing algorithm.
PACKMOL will then write the best solution to an MMP file.
The final product is special due to the distance between DGEBA and IPD molecules are now close enough to be cross-linked!!! Either another piece of vaporware could be written to cross-link the molecules or the code could be added to PACKMOL. The code will work as described in the "knitter" vaporware referenced in dynamic method one. The code will find the hot spots and the distances between molecules will be near perfect for cross-linking. The advantage of this approach is that LAMMPS does not need to execute to achieve cross-linking. The question is: is this really a dynamic method? Running LAMMPS in method one makes it dynamic. Running PACKMOL packing optimization is dynamic since multiple configurations of molecules are generated and cost functions calculated. Keep in mind that PACKMOL also allows for a pseudo-random seed to be specified allowing the researcher to generate multiple unique geometries for simulations. As mentioned in the "Cleaning House" article posted on 11 January, a structure cleaner would be a great piece of vaporware. I intend to test a structure cleaner against a final MD cell created using dynamic method two. This will be more complex given the number of molecules.
For those of you who read my blog please feel free to comment on my entries. I would love to hear your comments and please tell me where I am going wrong. I want to produce software that can be used. I welcome collaboration to create open source software for everyone. My next topic will address having to run LAMMPS prior to running LAMMPS. This is rather cryptic, but I have found that I must run LAMMPS on an initial MD cell to generate an initial MD cell at the correct density prior to going forward with simulations for scientific investigation.
structure monomer.mmp
number 1
center
fixed x,y,z,a,b,g
end structure
Obviously, either a new piece of vaporware is needed to generate the fixed monomer portions of a second PACKMOL input file or PACKMOL could be extended to make multiple passes in order to build up an MD cell.
In this case a second pass would recreate the monomers in fixed geometry and then the cross-linker molecules would be packed into spaces around them. An important question to ask about dynamic method one is where does the well mixed MD cell consisting of monomers and cross-linking molecules come from? As stated before the current PACKMOL packs the first molecule template followed by the second and so on. This segregates the molecules. If polymerization continues using dynamic method one we will eventually have a system that will not allow the diffusion of monomers or cross-linking molecules. Building a PACKMOL input file of alternating DGEBA and IPD molecules is doable with a Perl or Python script. An now back to our irregularly schedule blog entry...
PACKMOL uses constraints to put molecules as close together as possible without violating the user defined tolerance. 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. Remember, the molecular machine part file (MMP) created in NanoEngineer-1 can contain a "starting point", which I refer to as a "hot spot" If you think of a monomer in a polymer it has two hot spots. They are located at the sites where they react with other molecules. For my system of DGEBA and IPD these are at opposite ends of the monomer. DGEBA has two hot spots - male and female. IPD has four. IPD has two amine groups that have two hydrogens each. One IPD bonds with a maximum of four DGEBA. Therefore, if one wants to create a cross-linked epoxy system at 50 percent cross-linking one would specify only two of the IPD hot spots could be bonded and the number of IPD and DGEBA molecules would need to specified so that excess of either molecule does not occur - unless your system wants to study such a situation.
PACKMOL software changes could be moderate to complex. PACKMOL would pack molecules in a different manner than today. PACKMOL would have to understand the hot spots. When it encounters a hot spot it would use the coordinates of the hot spot to place the other molecule. The molecule to place would be based on the information stored in the hot spot when created in NanoEngineer-1. If there were a 50 percent probability of placing another monomer or a cross-linking molecule at the hot spot PACKMOL would flip a coin. If an atom has multiple hot spots and one has a better chance of reaction then PACKMOL will use this information stored in the hot spot. This new approach to accessing molecules for packing will not affect the packing algorithm. So what does this have to do with building cross-linked MD cells?
PACKMOL currently uses static constraints. One defines the MD cell dimensions and the tolerance between atoms of different molecules. One can select many 3-D MD cell shapes (e.g. spheres, boxes, etc...). Adding hot spots to PACKMOL will add a dynamic constraint. Every monomers would add two dynamic constraints as to where two IPD molecules have to exist. Specifying the placement of the hot spot of molecule B at tolerance distance from hot spot on molecule A sets up the necessary condition for cross-linking the two molecules. Keep in mind that PACKMOL will pack the molecules and maintain the value of the cost function over many attempts and then select the optimal packing configuration as the solution. Keep in mind that only the position of one hot spot is constrained while the attitude of the molecule to be place will vary in the packing algorithm.
PACKMOL will then write the best solution to an MMP file.
The final product is special due to the distance between DGEBA and IPD molecules are now close enough to be cross-linked!!! Either another piece of vaporware could be written to cross-link the molecules or the code could be added to PACKMOL. The code will work as described in the "knitter" vaporware referenced in dynamic method one. The code will find the hot spots and the distances between molecules will be near perfect for cross-linking. The advantage of this approach is that LAMMPS does not need to execute to achieve cross-linking. The question is: is this really a dynamic method? Running LAMMPS in method one makes it dynamic. Running PACKMOL packing optimization is dynamic since multiple configurations of molecules are generated and cost functions calculated. Keep in mind that PACKMOL also allows for a pseudo-random seed to be specified allowing the researcher to generate multiple unique geometries for simulations. As mentioned in the "Cleaning House" article posted on 11 January, a structure cleaner would be a great piece of vaporware. I intend to test a structure cleaner against a final MD cell created using dynamic method two. This will be more complex given the number of molecules.
For those of you who read my blog please feel free to comment on my entries. I would love to hear your comments and please tell me where I am going wrong. I want to produce software that can be used. I welcome collaboration to create open source software for everyone. My next topic will address having to run LAMMPS prior to running LAMMPS. This is rather cryptic, but I have found that I must run LAMMPS on an initial MD cell to generate an initial MD cell at the correct density prior to going forward with simulations for scientific investigation.
Friday, January 18, 2013
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:
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.
- 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.
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.
Friday, January 11, 2013
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.
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!
Thursday, January 10, 2013
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.
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.
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
Wednesday, January 2, 2013
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.
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.
Subscribe to:
Posts (Atom)