How to add two types of positive ions to a simulation box

[taoliu032]

Hi,

I wanted to add both NaCl and KCl to a simulation box by addSolvent(), but the method itself seems to only allow adding one type of positive ion at a time. I wonder if there is a clever way to work around this? I would guess this may be a common demand for others as well. I thought of two ways using PDBFixer.addSolvent(), but probably there is a better way. Please leave a comment if you have better idea(s)! Thanks!

Way1: add NaCl as normal first, then add KCl with 0nm as padding

fixer.addSolvent(padding=1.0*unit.nanometer, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0.14*unit.molar)
fixer.addSolvent(padding=0.0*unit.nanometer, positiveIon='K+', negativeIon='Cl-', ionicStrength=0.005*unit.molar)

The second step in this way will expand the initial water box a little bit: for example, box size goes from 64.5A to 65.9A and adds extra 4 water molecules. The expansion doesn't alter the [NaCl] too much, which is nice, but the extra space not being filled with water creates some kind of vacuum. When running simulation on it, the water molecules collapse a little and leave visible vacuum space in the box. I reckon it is not a good idea to leave vacuum in a simulation box, even with periodic boundary condition turned on.

Way2: add NaCl as normal, then remove water, then add KCl with 0nm as padding

fixer.addSolvent(padding=1.0*unit.nanometer, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0.14*unit.molar)
# remove the added water, only keep the NaCl
fixer.removeChains(3) # water chain is index 3
fixer.addSolvent(padding=0.0*unit.nanometer, positiveIon='K+', negativeIon='Cl-', ionicStrength=0.005*unit.molar)

In this way, when adding water in the second step, the resulting box water may be smaller than the initial box (the one that gets removed) -- because the NaCl may not spread across the whole initial box -- therefore can possibly strongly alter the effective [NaCl]. The bright side with this way is that the simulation box is filled with water (ie, added in the second step) and no vacuum will be created.

Any other better ideas?

[peastman]

I would just add all of it as NaCl, then transform half the Na atoms into K. Something like this should work (not actually tested, just to give the idea).

na = [r for r in topology.residues() if r.name == 'NA']
k = na[:len(na)//2]
for residue in k:
    residue.name = 'K'
    atom = next(residue.atoms())
    atom.name = 'K'
    atom.element = element.potassium

[taoliu032]

Oh that seems a lot easier! I was doing something similar, although by directly editing the PDB file. But pdbfixer removes some of the CONECT records, which is annoying. I didn't know how to edit the openmm topology directly -- now I know! I will give it a try, thanks Peter!

If it works, I will close the issue immediately :)

[taoliu032]

Yes, it works outrageously well! Thanks again!