One bad design decision I made was to use numerical atomtypes instead of the alphanumeric atomtypes the rest of the world uses. I did this to try to speed up achieving my masters thesis. At the time I had a lot of stuff happening in my life. Hopefully Summer of 2016 will be the Summer of Upgrades. The prototype changes for the WINDOWS uploaded to the SourceForge site.The new functionality is as follows:
1. Us the CFF91.frc file as provided with MSI2LMP as distributed with LAMMPS is now the standard.
2. PACKMOL now carries the alphanumeric atomtypes through from the MMP files to the final MMP file created.
3. The NanoEngineer-1 CAD tool now puts the alphanumeric atomtypes in the MMP file.
This upgrade includes the newest version of PACKMOL and MSI2LMP. I am working on building a
multithreaded PACKMOL EXE so we can make use of multi-threaded processors.
As always, if you use the tools please send me comments. I know I don't always answer. but I hope to solve that problem this summer.
I intend/hope to have Linux and OSX versions of PACKMOL and MSI2LMP upgraded by 15 July, 2016.
I have no time frame for a Linux NanoEngineer-1 yet. This has to be built from a different source tree and suffers more from and tool ignorance.
Sunday, July 3, 2016
Tuesday, December 8, 2015
Constructing Systems of Molecules for Molecular Dynamics Simulations Using Self Assembly
It's December 2015 and I have spent many an hour reading about how others build systems of molecules for molecular dynamics simulations. There are many ways to Rome and all have obstacles and suffer from problems and seem not to provide good, computationally quick and cheap solutions. My biggest problem since last October has been funding. I solved this problem in February of 2015 by going back to work full-time. Sadly, working full-time leaves much less time to devote to my research. This is why I have spent time reading papers on building systems. I have seen 12 step programs and and 21 step programs to build a system. The definition of a system is one or more molecules in equilibrium prior to running molecular dynamics (MD) to determine some physical property or perhaps to collect data for diffusion of small molecules in a polymer material. Digging deeper the monomers must be connected to form the polymer. Polymer chain lengths must be decided and constructed and the bonds between polymers must be constructed. The bonds between polymer chains are cross-links. This is where the fun begins. The number of and types of steps seem predicated on achieving some desired outcome for the initial simulation system. Outside of the theoretical (e.g. first principles) work of Flory (Rotational Isomeric State Theory), Theodorou and Suter ("Detailed Molecular Structure of a Vinyl Polymer Glass) the phenomenological concepts continue to rule system construction and therefore "free volume" concepts are invoked.
Constructing a system is usually defined as the automation of the polymerization process and/or constructing bonds between polymers, where molecules are attached to one another. In the case of an epoxy system, DGEBA/IPD, one IPD molecule can cross-link with up to DGEBA molecules. This can result in a 3D networked system of monomers. Many of the reviewed approaches literally beat the system to death in terms of cycling temperature and pressure of the system based on chemistry of the system until it exhibits experimental properties of the macroscale system in simulation. The usual suspects compare the built system to experimental data via radial distribution functions and the like and declare success. How can one compare a simulated system with an experimental system? I suggest comparing elastic, shear, bulk moduli and density values to their experimental counterparts at specific temperature ranges. Using molecular dynamics (MD) one must average over some significant time duration to give an average (time independent) value. One can do something similar with Monte Carlo methods (MCM), but it does not contain the kinetic energy (KE). Based on many sources across many disciplines, one relies on the lowest system energy representing the most probable configuration, where configuration is geometry.
Obviously, geometry and potential energy (PE) are related via force fields. MD subjects the simulated system to both potential and kinetic energy (KE) since MD provides velocity. Experimental physical properties such as moduli, density and glass transition temperature are are based on PE and KE and therefore dynamic. The interesting point is PE is usually negative and KE is positive. Looking at the equations for the CFF force fields like COMPASS one sees that a single molecule can have a zero PE value if every bond length, bond angle, out-of-plan and dihedral angles are all at their equilibrium values.Once energy is stored by distorting bond lengths and/or angles the system energy is decreased. The movement of atoms provide the kinetic energy, which increases the system energy. Of course multiple molecules create a system that must account for Van der Waals and Coulombic forces as COMPASS does. Forces between atoms are time varying and atoms are attracted and repelled based on proximity. If one were to take a snapshot in time of a system of molecules one would see atoms in all molecules at positions where bond lengths and angles away from the zero PE positions. Assuming that none of the molecules have a net delta velocity the system can be considered to be in equilibrium. Considering a mass of epoxy in an ocean of nitrogen molecules exerting 1 bar of pressure the epoxy would consist of bond lengths and angles dictated by the external forces. However, atoms in a material at some temperature must move as per normal modes of vibration with stretching, bending, rocking, and wagging motions. The result is no net translation or rotation. For molecules in close proximity their Van der Waals and Coulombic forces create a viscous-like environment to dampen movement.
Obviously, in the real-world cross-links form in a nonequilibrium environment with mass transfer and most likely temperature gradients. Chemical reaction can be endo and/or exothermic. Molecular dynamics doesn't account for electronic configurations. There are some reactive force fields, but they are specific to certain types of chemical reactions. Barring advances in this area the question remains how does one get molecules to "assemble" themselves into the lowest energy conformations? So far, the literature cites using a proximity strategy. Once two molecules come close enough (e.g. 6-10 Angstroms) at their reactive sites a bond is created between the two molecules. This is the cross-link bond. The two molecules are now one molecule and the PE of the one molecule depends on the rotation of the smaller molecule around the newly formed bond and the length of the newly formed bond assuming the two original molecules don't collide. A human can connect the two molecules and with a little playing "eyeball" where the two molecules should be oriented to minimize the PE. Running a MMEM software will attempt to minimize the energy of the entire system. If the two molecules are bonded and the new bond forms deviates from the zero PE bond length or if the angles it forms deviate from the zero PE angles then the new molecule increases in PE. Time must be spent relaxing the system to a lower energy and there is no guarantee of moving the system to the best conformation. In summary, molecules need to move relative to one another, when they are in proximity they must be bonded together forming the cross-link and then the system PE must be adjusted. Proximity is a bad strategy for cross-linking.
How do molecules move? The molecules or the atoms composing molecules can move by dynamic, static and energy methods. If one adds proximity cross-linking to these molecule movement methods one defines possible system construction methods. The current construction methods as follows:
1. Dynamic Methods - Molecular Dynamics (MD) is obviously a time series of atom and molecule
movements in a confined space. Molecules go whizzing around in space. force fields (think
springs between atoms) are used to insure the molecules do actually translate and rotate according to Newton's laws. Proximity is used to cross-link molecules at reaction sites to form cross-links.
2. Static Methods - Self-avoiding walks (SAW) and Packmol optimization methods fall into
this category. Self-avoiding walks place atoms and test for collisions and backtrack when
necessary and then try another direction to place atoms. Packmol calculates a packing
optimization solution based on a penalty function for atoms with respect to one another and to
geometric boundaries (e.g. walls ). Molecules are moved via a heuristics algorithm. Proximity is
used to cross-link molecules at reaction sites to form cross-links.
3. Energy Methods - Metropolis Monte Carlo Method (MMCM) and Simulated Annealing (SA) fall
into this category. SA relies on reducing the system temperature of the system over the course of the computation. The probability of the small spatial move of a molecule is calculated accounting
for the change in system energy and compared to a random number. Metropolis Monte Carlo
Method is done with a single system temperature. The methods are not time-based and in many
case higher energy moves are accepted that allow the system to bypass high energy barriers and
hopefully find a lower system energy. Proximity is used to cross-link molecules at reaction sites
to form cross-links.
4. Other Method(s) - Genetic Algorithms, reactive force fields, etc...
Packmol uses a penalty function to determine if the packing of molecules is successful. Packmol randomly distributes molecules in a given volume and then employs a linked-cell approach to determine if every atom of a molecule is at least 2 Angstroms from the atoms of all other molecules. The linked-cell approach linearizes this problem since only atoms in the linked-cell need to be compared drastically reducing the number of comparisons. Packmol uses heuristics to move molecules with large penalty functions to locations where molecules with low penalty functions are located. The linked-cell and proximity is just part of the penalty function calculation. Every atom has a list of restrictions. A simple system contains a single box restriction where all atoms must be inside the box. Unless the atom is very near the perimeter its restriction contribution is zero. An atom alone in a linked-cell has the same penalty function value. Simple restrictions such as the box boundary is read in from a configuration file. Every atom has a component contribution of the penalty function. Molecules with the highest penalty function are too close to other molecules and are translated and rotated close to molecules with the lowest penalty function values.
GENCAN is used to minimize the penalty function described above. This is done by calculating the values based on atom proximity and restrictions. The penalty function is a function of the translation and rotation of the molecules in the system. A penalty function of zero guarantees no spearing, or bad dynamics for the system of molecules.
When the atom proximity is being calculated atom types can be used to match the correct pair of atoms across two molecules. In the case of the DGEBA/IPD system the carbon and nitrogen atoms need to be some distance apart to support cross-linking (6-10 Angstroms). If this is the only criteria then bad geometry occurs due to the two molecules involved are not constrained enough to force GENCAN to consider only the dihedral angle between the two molecules. Also, the currently implemented heuristics often forces the two molecules in question apart making the cross-link much greater in length forcing the cross-link to be abandoned. Ideally, two molecules (e.g. IPD and DGEBA) are in proximity and a new restriction could be created that gets calculated in the restrictions section of the algorithm. This on-the-fly restriction (or dynamic constraint) now states that two atoms in different molecules must be some distance apart. As molecules move now the molecule that is to be cross-linked to it must move as well. The solution to this problem comes from the geometry of the sp3 nitrogen and the sp3 carbon involved in the cross-link site between the two molecules. The nitrogen has two hydrogens which come off during the reaction with the epoxide ring which provides the sp3 carbon. The key is if the carbon were to be located where one of the hydrogens are and if the nitrogen were to be located where one of the hydrogens is on the sp3 carbon then the nitrogen and carbon form an axis and the two molecules have the correct attitude relative to one another with the exception of their dihedral angle. By not fixing the dihedral angle, the molecules are free to rotate about their common axis using the GENCAN algorithm. This new approach developed is to specify the hydrogens that are not needed after cross-linking as having an atom type of zero. This is jokingly referred to as "zeronium." The "zeronium atoms are later removed when the molecular new system's machine part file is written. While they exist in the packing process they provide a coordinate to locate the other molecule's Carbon or Nitrogen atoms. The goal being to achieve the correct geometry in terms of bond lengths and bond angles initially and letting GENCAN solve the dihedral angle. Linking together molecules in this manner solves the problem of creating new molecules and allowing the cross-linked molecules to achieve new conformations over the course of the packing optimization. It solves the bad geometry problem by taking the ideal geometry as a first estimate.
The heuristics algorithm become very different. The new algorithm becomes a process of "self assembly." If each molecule were a vehicle capable of navigating in 3D they could broadcast a need for other molecules to attach themselves to them. The attachment site coordinates are the data needed for an prospective molecule to change it's attitude in space to line up for attachment. This approach leads to a self assembly heuristics algorithm. Since two points in space are used to cross-link two molecules, the GENCAN algorithm is free to solve the dihedral angle. With the dynamic constraints moving a single molecule would require the cascading movement of multiple molecules. This is why moving free monomer molecules takes priority. Proximity is tested for during the GENCAN algorithm and dynamic constraints are created on the fly as well.
The current Packmol software has no concept of "holes" or "voids" inside the box. Also, one successfully runs the software on successively smaller and smaller boxes until the algorithm fails to pack the molecules. Making the box bigger gets one back to a successful packing. This usually provides the smallest possible box with the least amount of "holes" in it. The new heuristics will require moving molecules to a new location in the box with the purpose of creating cross-links. The cross-link sites on each molecule are known and moving the IPD into proximity of the DGEBA or the otherway around will in fact create the permanent linkage between the molecules while the dynamic constraint that follows provides the correct geometry.
The modifications to be implemented use the proximity penalty function algorithm to identify cross-link pairs. The cross-links are maintained in a new dynamic data structure. At the same time the new dynamic constraint is created for the carbon atom in the DGEBA molecule and the zeronium atom in the IPD molecule. The zeronium atom in the DGEBA is paired up with the nitrogen in the IPD molecule. The number of restrictions per atom is currently set at 10 and can be increased to handle these dynamic constraints. Unlike MD software, Packmol doesn't care when two atoms occupy the same space in this case it helps define the necessary geometry and to avoid overlaps and collisions. The coding of the Packmol enhancement is in progress and shows promise in self assembling the monomers to polymers and in creating cross-links.
Constructing a system is usually defined as the automation of the polymerization process and/or constructing bonds between polymers, where molecules are attached to one another. In the case of an epoxy system, DGEBA/IPD, one IPD molecule can cross-link with up to DGEBA molecules. This can result in a 3D networked system of monomers. Many of the reviewed approaches literally beat the system to death in terms of cycling temperature and pressure of the system based on chemistry of the system until it exhibits experimental properties of the macroscale system in simulation. The usual suspects compare the built system to experimental data via radial distribution functions and the like and declare success. How can one compare a simulated system with an experimental system? I suggest comparing elastic, shear, bulk moduli and density values to their experimental counterparts at specific temperature ranges. Using molecular dynamics (MD) one must average over some significant time duration to give an average (time independent) value. One can do something similar with Monte Carlo methods (MCM), but it does not contain the kinetic energy (KE). Based on many sources across many disciplines, one relies on the lowest system energy representing the most probable configuration, where configuration is geometry.
Obviously, geometry and potential energy (PE) are related via force fields. MD subjects the simulated system to both potential and kinetic energy (KE) since MD provides velocity. Experimental physical properties such as moduli, density and glass transition temperature are are based on PE and KE and therefore dynamic. The interesting point is PE is usually negative and KE is positive. Looking at the equations for the CFF force fields like COMPASS one sees that a single molecule can have a zero PE value if every bond length, bond angle, out-of-plan and dihedral angles are all at their equilibrium values.Once energy is stored by distorting bond lengths and/or angles the system energy is decreased. The movement of atoms provide the kinetic energy, which increases the system energy. Of course multiple molecules create a system that must account for Van der Waals and Coulombic forces as COMPASS does. Forces between atoms are time varying and atoms are attracted and repelled based on proximity. If one were to take a snapshot in time of a system of molecules one would see atoms in all molecules at positions where bond lengths and angles away from the zero PE positions. Assuming that none of the molecules have a net delta velocity the system can be considered to be in equilibrium. Considering a mass of epoxy in an ocean of nitrogen molecules exerting 1 bar of pressure the epoxy would consist of bond lengths and angles dictated by the external forces. However, atoms in a material at some temperature must move as per normal modes of vibration with stretching, bending, rocking, and wagging motions. The result is no net translation or rotation. For molecules in close proximity their Van der Waals and Coulombic forces create a viscous-like environment to dampen movement.
Obviously, in the real-world cross-links form in a nonequilibrium environment with mass transfer and most likely temperature gradients. Chemical reaction can be endo and/or exothermic. Molecular dynamics doesn't account for electronic configurations. There are some reactive force fields, but they are specific to certain types of chemical reactions. Barring advances in this area the question remains how does one get molecules to "assemble" themselves into the lowest energy conformations? So far, the literature cites using a proximity strategy. Once two molecules come close enough (e.g. 6-10 Angstroms) at their reactive sites a bond is created between the two molecules. This is the cross-link bond. The two molecules are now one molecule and the PE of the one molecule depends on the rotation of the smaller molecule around the newly formed bond and the length of the newly formed bond assuming the two original molecules don't collide. A human can connect the two molecules and with a little playing "eyeball" where the two molecules should be oriented to minimize the PE. Running a MMEM software will attempt to minimize the energy of the entire system. If the two molecules are bonded and the new bond forms deviates from the zero PE bond length or if the angles it forms deviate from the zero PE angles then the new molecule increases in PE. Time must be spent relaxing the system to a lower energy and there is no guarantee of moving the system to the best conformation. In summary, molecules need to move relative to one another, when they are in proximity they must be bonded together forming the cross-link and then the system PE must be adjusted. Proximity is a bad strategy for cross-linking.
How do molecules move? The molecules or the atoms composing molecules can move by dynamic, static and energy methods. If one adds proximity cross-linking to these molecule movement methods one defines possible system construction methods. The current construction methods as follows:
1. Dynamic Methods - Molecular Dynamics (MD) is obviously a time series of atom and molecule
movements in a confined space. Molecules go whizzing around in space. force fields (think
springs between atoms) are used to insure the molecules do actually translate and rotate according to Newton's laws. Proximity is used to cross-link molecules at reaction sites to form cross-links.
2. Static Methods - Self-avoiding walks (SAW) and Packmol optimization methods fall into
this category. Self-avoiding walks place atoms and test for collisions and backtrack when
necessary and then try another direction to place atoms. Packmol calculates a packing
optimization solution based on a penalty function for atoms with respect to one another and to
geometric boundaries (e.g. walls ). Molecules are moved via a heuristics algorithm. Proximity is
used to cross-link molecules at reaction sites to form cross-links.
3. Energy Methods - Metropolis Monte Carlo Method (MMCM) and Simulated Annealing (SA) fall
into this category. SA relies on reducing the system temperature of the system over the course of the computation. The probability of the small spatial move of a molecule is calculated accounting
for the change in system energy and compared to a random number. Metropolis Monte Carlo
Method is done with a single system temperature. The methods are not time-based and in many
case higher energy moves are accepted that allow the system to bypass high energy barriers and
hopefully find a lower system energy. Proximity is used to cross-link molecules at reaction sites
to form cross-links.
4. Other Method(s) - Genetic Algorithms, reactive force fields, etc...
Packmol uses a penalty function to determine if the packing of molecules is successful. Packmol randomly distributes molecules in a given volume and then employs a linked-cell approach to determine if every atom of a molecule is at least 2 Angstroms from the atoms of all other molecules. The linked-cell approach linearizes this problem since only atoms in the linked-cell need to be compared drastically reducing the number of comparisons. Packmol uses heuristics to move molecules with large penalty functions to locations where molecules with low penalty functions are located. The linked-cell and proximity is just part of the penalty function calculation. Every atom has a list of restrictions. A simple system contains a single box restriction where all atoms must be inside the box. Unless the atom is very near the perimeter its restriction contribution is zero. An atom alone in a linked-cell has the same penalty function value. Simple restrictions such as the box boundary is read in from a configuration file. Every atom has a component contribution of the penalty function. Molecules with the highest penalty function are too close to other molecules and are translated and rotated close to molecules with the lowest penalty function values.
GENCAN is used to minimize the penalty function described above. This is done by calculating the values based on atom proximity and restrictions. The penalty function is a function of the translation and rotation of the molecules in the system. A penalty function of zero guarantees no spearing, or bad dynamics for the system of molecules.
When the atom proximity is being calculated atom types can be used to match the correct pair of atoms across two molecules. In the case of the DGEBA/IPD system the carbon and nitrogen atoms need to be some distance apart to support cross-linking (6-10 Angstroms). If this is the only criteria then bad geometry occurs due to the two molecules involved are not constrained enough to force GENCAN to consider only the dihedral angle between the two molecules. Also, the currently implemented heuristics often forces the two molecules in question apart making the cross-link much greater in length forcing the cross-link to be abandoned. Ideally, two molecules (e.g. IPD and DGEBA) are in proximity and a new restriction could be created that gets calculated in the restrictions section of the algorithm. This on-the-fly restriction (or dynamic constraint) now states that two atoms in different molecules must be some distance apart. As molecules move now the molecule that is to be cross-linked to it must move as well. The solution to this problem comes from the geometry of the sp3 nitrogen and the sp3 carbon involved in the cross-link site between the two molecules. The nitrogen has two hydrogens which come off during the reaction with the epoxide ring which provides the sp3 carbon. The key is if the carbon were to be located where one of the hydrogens are and if the nitrogen were to be located where one of the hydrogens is on the sp3 carbon then the nitrogen and carbon form an axis and the two molecules have the correct attitude relative to one another with the exception of their dihedral angle. By not fixing the dihedral angle, the molecules are free to rotate about their common axis using the GENCAN algorithm. This new approach developed is to specify the hydrogens that are not needed after cross-linking as having an atom type of zero. This is jokingly referred to as "zeronium." The "zeronium atoms are later removed when the molecular new system's machine part file is written. While they exist in the packing process they provide a coordinate to locate the other molecule's Carbon or Nitrogen atoms. The goal being to achieve the correct geometry in terms of bond lengths and bond angles initially and letting GENCAN solve the dihedral angle. Linking together molecules in this manner solves the problem of creating new molecules and allowing the cross-linked molecules to achieve new conformations over the course of the packing optimization. It solves the bad geometry problem by taking the ideal geometry as a first estimate.
The heuristics algorithm become very different. The new algorithm becomes a process of "self assembly." If each molecule were a vehicle capable of navigating in 3D they could broadcast a need for other molecules to attach themselves to them. The attachment site coordinates are the data needed for an prospective molecule to change it's attitude in space to line up for attachment. This approach leads to a self assembly heuristics algorithm. Since two points in space are used to cross-link two molecules, the GENCAN algorithm is free to solve the dihedral angle. With the dynamic constraints moving a single molecule would require the cascading movement of multiple molecules. This is why moving free monomer molecules takes priority. Proximity is tested for during the GENCAN algorithm and dynamic constraints are created on the fly as well.
The current Packmol software has no concept of "holes" or "voids" inside the box. Also, one successfully runs the software on successively smaller and smaller boxes until the algorithm fails to pack the molecules. Making the box bigger gets one back to a successful packing. This usually provides the smallest possible box with the least amount of "holes" in it. The new heuristics will require moving molecules to a new location in the box with the purpose of creating cross-links. The cross-link sites on each molecule are known and moving the IPD into proximity of the DGEBA or the otherway around will in fact create the permanent linkage between the molecules while the dynamic constraint that follows provides the correct geometry.
The modifications to be implemented use the proximity penalty function algorithm to identify cross-link pairs. The cross-links are maintained in a new dynamic data structure. At the same time the new dynamic constraint is created for the carbon atom in the DGEBA molecule and the zeronium atom in the IPD molecule. The zeronium atom in the DGEBA is paired up with the nitrogen in the IPD molecule. The number of restrictions per atom is currently set at 10 and can be increased to handle these dynamic constraints. Unlike MD software, Packmol doesn't care when two atoms occupy the same space in this case it helps define the necessary geometry and to avoid overlaps and collisions. The coding of the Packmol enhancement is in progress and shows promise in self assembling the monomers to polymers and in creating cross-links.
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.
Subscribe to:
Posts (Atom)