[gmx-users] Some Issues on Substrates/PBC

[Michael Hirtz]

Dear Listmembers,

This post is quite lengthy because it contains several different problems that
I came across during my simulations and that I cannot figure out up to now. I
included some links to pictures to give you some visualization and would be
glad for every suggestion or hint even to part of my problems.

Generally for my simulations I need substrates of silicon oxide and graphite
(each one in different simulations). The substrate fills the simulation box
completely in x-y direction so it should look like an infinite surface
considering the periodic boundary conditions.

The simulation with silicon oxide is working fine but I experienced some
strange thing though: When the .gro file contains negative coordinates and I
generate a .tpr for 4 nodes the structure is messed up totally after energy
minimization (although it converges and gives no errors). I guess that the
breaking up in four parts is somehow connected to the 4 nodes, because
generating and running a 1 node .tpr solves this problem and the structure
keeps intact after em - see pictures here: ( http://tinyurl.com/59zfda ). When
I shift the coordinates in the .gro with editconf so that every coordinate is
positive the structure can be em with both 1 or 4 node without problems.

For the main md I have to stick with 1 node anyway because I get a "Shake
block crossing node boundaries" error when I try to use 4 nodes. I found the
advise to use less nodes in the mailing list archives and as I said it works
for me, too, but is there any other way to get around this? For the test runs
I'm doing now 1 node is fine, but for the real simulations later on I would
like to use more nodes because the system will get bigger...

With the graphite I have another problem: Again I have to stick with only one
node, otherwise the structure is broken after energy minimization (with 4
nodes http://tinyurl.com/5eqngz ). Additionally here the energy minimization
is not converging even after 2500 steps, forces keep in the range of 10^4 but
I can still do a md with this file that is not exploding (md with the starting
structure without any energy minimization will result in an exploding system
as one might expect). After energy minimization the graphite sheets (four are
in the simulation) are bended up/down correspondently to the connected
periodic boundary at the edges of the simulation box (box is triclinic,
http://tinyurl.com/5m72ke ). When I run an md with this as starting structure
without position restraints on the carbon atoms, the sheets cleave beginning
from the edges and then reunite in a curled conformation within 10ps (
http://tinyurl.com/5j4l3z ).

If I put position restraints on all carbon atoms the graphite keeps flat, but
of course I still have the up/down bending at the box edges, which I also want
to get rid off to have a real "infinite" graphite surface.

First I thought that this was maybe a problem of my graphite structure not
fitting into my simulation box (hence the curling of the surface) but I
checked the bonding lengths and b0 in the topology, the unit cell vectors and
box vectors thoroughly and everything seems to be right. If I copy my
simulation box and shift it by one of the box vectors I get a perfect
continuation without any deviation on the boundaries, so everything seems like
it should be for the starting structure.

I also simulated the same graphite sheets within a bigger box, so that border
carbons could not feel their counterparts over the periodic boundaries. That
gave no distortion at the edges after em and the em converged, but of course
in this way I cannot get an continous surface.

Even stranger is, that after em the whole structure coordinates are shifted
into the neighbouring periodic cell! For a bigger system with some
additionally molecules that should absorb to the graphite surface the graphite
sheets even ended up in different periodic images of the simulation box after
em ( http://tinyurl.com/5uaxfa ). Maybe that's more a visualization problem
because the following md simulation shows, that the molecules adsorb to the
surface nevertheless (as one would suggest from pbc). Of course I can
reconfigure my .gro for visualization by substracting or adding the box
vectors untill all atoms end up inside the 'real' simulation box, but is it
really supposed to behave that way, or is there something definitely wrong
with my structure?

I saw this "shifting" of the whole substrate into another periodic image of
the simulation box also in the simulation with the silicon substrate: it
'jumped' into the periodic cell and then jumped back one frame later.
The other present (small) molecules stayed in the real simulation box all the
time. I guess that this is maybe a problem caused by the substrate "molecule"
filling up the whole simulation cell (at least in xy), but is there a way of
getting Gromacs to always write out coordinates that lie within the real
simulation cell instead of an periodic image?

And is there a way of defining bonds over the periodic boundaries? When I
tried this it always seemed as if Gromacs interpreted these bonds as extremly
overstretched bonds (that are stretched over the whole simulation box
instead of just reaching the "short way" over the periodic boundary) yielding
crippled structures ( http://tinyurl.com/54uvs4 ).

All simulations used the ffG43b1 force field. Cutoffs of 1.2 and an Rlist of
1.4. Coulomb was shifted from 0.0 and VdW from 0.9. Timestep was 2 fs and pbc
= xyz. The boxes were large in z direction to generate xy pbc. I can also put
the topologies online if someone would like to look at them.

I hope anyone can give me some advice because although I tried a lot of things
right now I feel not like coming  close to a solution for my problems. :(

Best regards,
    Michael
-- 
http://www.defux.de