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.