I want to insert multi-atom anions such as sulphate (SO4), carbonate (CO3), and thiocyanate (SCN-) into an existing lammpsdata file. The lammpsdata file contains a single oligomer along the body diagonal of a cube and is surrounded by water molecules. I have built this using moltemplate using the OPLSAA forcefield parameters. Unfortunately, moltemplate doesn’t have the OPLSAA parameters for these ions so I cannot add them that way. I have prepared xyz, pdb and lammpsdata files for these ions. From here I have multiple ideas on how to proceed with the insertion, but I don’t know which would be the proper way:
Read the lammpsdata files for the ions after the main read_data command. The problem is controlling the location of the insertion. I think it’s an extremely tedious job to proceed this way.
Instead of creating the simulation box directly in moltemplate, create only the polymer lammpsdata file so that moltemplate gives me the OPLSAA files for LAMMPS (pair bond dihedral and angle coefficeints, charges etc.). Now I would have to convert this lammpsdatafile to xyz format and assemble my polymer, ions, and water molecules in PACKMOL.
I am thinking the 2nd method is more sensible, but I am still unsure how to use packmol correctly. Could anyone please guide me into preparing my simulation cell correctly using PACKMOl. Also, I’m afraid that PACKMOL will change the atom type definition assigned to my polymer by moltemplate. I have been reading many research papers but not even a single one describes how to build the inital system.
I am also attaching the lammpsdata file that I have built using moltemplate.
I disagree. The problem is that by converting your complete topology back to an .xyz or .pdb file you lose all information about topology and atom type assignment.
Your problem is certainly a difficult one.
I see two additional and more promising ways to approach this:
create suitable molecule files for your ions and then use fix gcmc to insert them (you may need to use soft-core potentials to avoid divergence of forces when MC attempts are made where atoms come to close which would lead to invalid math operations). The problem here is that you want to insert charged objects, but this kind of procedure would work best for neutral objects, i.e. you would need to also insert counter ions. It would be worse, if the neutralizing charges are already in the system. One way to address this could be to temporarily neutralize the existing system and the anions until you have the desired composition (it is easy to remove one or more excess anions later to achieve that state) and only then restore the original charge assignments and continue equilibrating. Since these steps are before doing a proper equilibration, you can take a lot of liberties in terms of using settings that would be “wrong”, but will accelerate the process of reaching the system composition you desire.
temporarily change your pair style to “soft” and ignore that there are overlaps and just add the anions anyway, but do not compute regular OPLS forces at first. Instead, you should use fix adapt to start with no non-bonded interactions at first and then gradually “grow” those over the course of the simulation until the conflicts due to overlaps have been resolved. Please see the “micelle” example in the LAMMPS source code distribution for a (very simple) case where this is done.
In general, any type assignment (atom, bond, angle, dihedral, etc.) and bond topology information is valuable information and thus must be maintained, as it will be very difficult to recover it after it has been removed. It is easier to modify the interactions so they don’t cause problems.
Thanks for the clarifications and suggestions for proceeding with the insertion. I still have a few questions regarding these suggestions. I am trying to implement this method:
temporarily change your pair style to “soft” and ignore that there are overlaps and just add the anions anyway, but do not compute regular OPLS forces at first.
When you say this do you mean that I need to use fix gcmc to add the anions anyway?
Also, when you say that I shouldn’t compute any OPLS forces does it mean that I need to rewrite all my pair_coeff temporarily?
If I were to implement fix gcmc for insertion, I have been trying to follow the example given in examples/mc but I am unable to replicate it for my system. I have never used fix gcmc before so I might be doing it wrong.
I am attaching the files that I have used for my attempt PhnRvTjN-Anions.zip (1.2 MB)
In this the run.*.lmp is my central input script file.
Here it looks like the user Amalie was adding the molecules simply to a region by using the create_atoms command.
So, I am wondering if I could purposefully leave out voids while creating my lammpsdata file using moltemplate. I can mark these void regions using the region command and then use the create_atoms command as implemented in this thread. I know this can be a tedious attempt but it certainly feels less intimidating than having to mess up my pair_coeff definitions.
Could you please shed more light on my approach. I seem to have a hard time trying to get this thing to work.
You can use write_data ... nocoeff and write_coeff to separate topology and force field parameters. That will make it easy to switch back the pair style and restore the existing coefficient settings via include.
Sorry, but I don’t have the time to review, debug, and correct your inputs. If anything that would be the job of your adviser/tutor. Basically, I can point you in the direction, but you have to figure out the rest.
But it has conceptual problems. You want random fairly random placement of your anions. If anything, you want to find locations for your anions that have a favorable electrostatic potential. VMD has a cionize tool that can do that for simple (atomic) ions, which can be a bit time saver for time to equilibrium, if you have to place ions to neutralize a DNA strand or similar.
You are tackling a “hard” problem. Some complications are to be expected.
My preferred option would still be the switch to a soft (core) potential and then gradually grow the potential so that atoms will start repelling each other and thus gently unoverlap them over time. If you need to preserve the secondary structure for some part of the system, you can always define a group with these atoms and apply fix spring/self to that.
TL;DR, there are many tools in LAMMPS that can be useful for your case, but it requires some thinking and experimenting and the guidance of an experienced adviser to navigate through the potholes without too much damage. As with most simulations worth doing, the major effort is not so much running the simulation, but the planning, the setup and the analysis.
Prepare your “ingredients list”, preferably in some kind of log book you can refer back to. Order matters. For example:
1 oligomer
4004 H2O
50 Na
25 SO4
Prepare all single molecules as pdb or xyz inputs for PACKMOL. I’d recommend stripping the waters from your existing data file and letting PACKMOL place them as well, but you could also keep them together as is. Also prep single molecule data files – mainly because they will be useful for building the topology in step 6.
Prepare your list of atom types – that is, types 1 to 18 will be in the oligomer, 19 will be a water hydrogen, 20 will be a water oxygen, etc. Preparing them in the order that molecules are listed in step 1 may help you avoid some errors.
Build your box in PACKMOL, putting your molecules into the input in the order you decided in step 1. Get your output as xyz.
Build your Atoms section with columns for atom ID, molecule ID, atom type, and charge, and the last three columns for XYZ coordinates from the PACKMOL output. This can be relatively straightforward to build in awk. You can also use Bash for loops to build the first four columns and paste to merge it with the PACKMOL coordinates.
Build your Bonds section, reusing your existing data file for the oligomer and programmatically generating the rest. Repeat for angles, dihedrals, and impropers. Your list from steps 1 and 2 will be invaluable here.
Merge together all your generated files. You can reuse all the prior sections from the oligomer data file, updating all the numbers and numbers of types as needed.
Is this all tedious? Yes. You can probably save yourself a lot of time and error by using Moltemplate to generate the topology (step 6), and it will accept XYZ input for setting coordinates, making step 5 much easier too. But I’ve written out this procedure to tell you how you could get your desired result using only PACKMOL and basic text processing tools.
I would follow the procedure described by Shern. Since you need to define the topology and parameters of the extra ions anyway, it makes sense to me to create a LAMMPS template file with the definition of each ion type and use Moltemplate to render the initial guess (XYZ) plus updated force field.
You can use the existing OPLSAA template to learn the syntax and create a new file with the extra molecules. It is fairly easy to create the topology of a single molecule, e.g. from a Z-matrix (use Molden to create one).