[gmx-users] Deprotonation and charging of a non-protein in salt solution

[Mark Abraham]

Steven Kirk wrote:
> Hello,
> 
> I have a large (5nm diameter) particle of amorphous silica in my 
> simulation box, generated using x2top and some OPLS-AA parameters and 
> kind advice from David Spoel. The particle consists of only Si, O, and 
> H, and all atoms show up in my topology as uncharged. All of the H atoms 
> are hanging off Si-OH groups on the particle surface.
> 
> I now want to partially deprotonate the particle, to mimic the effect of 
> pH in a TIP4P water solvent box with a known fixed number N of randomly 
> inserted Na+ counterions present(inserted using genion), so that the 
> particle acquires a negative charge and the simulation cell is 
> electrically neutral.
> 
> My naive strategy looks like this:
> 
> 1. Randomly select a number N of H atoms on the particle surface and 
> remove them.
> 2. Set a charge of -1 on those surface O atoms which have been 
> deprotonated (how?).
> 3. Let the whole cell contents relax (eventually, hopefully, to an 
> particle surrounded by a stable electrical double layer)
> 
> Can anyone suggest an easy (preferably automated) procedure to generate 
> such a (hopefully chemically consistent) configuration? What files would 
> need to be edited and how? Even just one such possible configuration 
> would be useful, but of course the ability to generate an ensemble of 
> such randomly edited structures progammatically would be even more useful.

I'd expect you will have to write some code to do it.

Given a pre-existing silica particle topology, I would define a virtual 
sphere centered on the same point and choose three deprotonation points 
randomly on it (easier said than done) such that they were not all on 
the same great circle. Now, in principle, you can choose further 
deprotonation points on the surface of the sphere to maximise the 
all-against-all distance sum between deprotonated points (through-space 
distance, not sphere surface distance). You would need some kind of 
optimization algorithm to seek that maximum, and have to try a few 
random starting points to make sure you get the global optimum at each 
step. This means when you've made N deprotonation points, they aren't 
clustered together, which would be unphysical. (genion probably has a 
suitable algorithm, if you swap signs of charges and can constrain it to 
a sphere...)

Now you need a (for example, Perl or Python) script to map the 
deprotonation points on the sphere to the nearest hydroxides and convert 
them to oxides and adjust the topology files accordingly. Obviously the 
oxygen will change its atom type.

Now you've got a particle that's been deprotonated in a hand-wavingly 
reasonable way, you can use genion to plunk down enough Na+ outside the 
particle to neutralise, and away you go.

Mark

P.S. Give it to an undergrad for a project :)