GROMACS version: 2021
GROMACS modification: No
Here post your question
Dear all,
I am performing a Computational Electrophsiology simulation between 2 membrane. I want to see how the potassium ion crossing the channel and using the swapping to ensure A side has more K ion than the B side. However, this error come up :
Fatal error:
Could not get index of W atom. Compartment contains 152 W molecules before
swaps.
However if I set the opposite, i.e. B side has more ion than A side, the system run a little bit before blowing up… this is the error:
Warning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.
Box (3x3):
Box[ 0]={Warning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.
-nan}
Box[ 2]={ -nan, -nan, -nan}
Can not fix pbc.
Box (3x3):
Box[ 0]={ -nan, , -nan 0.00000e+00, , 0.00000e+00 Box[ 0]={ Box[ 0]={, -nan , Box[ 0]={ 0.00000e+00 -nan, Warning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.
Box (3x3):
Box[ 0]={ -nanWarning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.
hence I am a bit lost what is happening and how do I solve this?
Below is how I set the .mdp up:
swapcoords = Z ; Swap positions: no, X, Y, Z
swap-frequency = 50 ; Swap attempt frequency
split-group0 = layer1 ; Defines compartment boundary
split-group1 = layer2 ; Defines other compartment boundary
massw-split0 = no ; use mass-weighted center?
massw-split1 = no
solvent-group = W ; Group containing the solvent molecules
iontypes = 2 ; Number of different ion types to control
iontype0-name = K
iontype0-in-A = 1751 ;
iontype0-in-B = 64
iontype1-name = CL
iontype1-in-A = -1
iontype1-in-B = -1
coupl-steps = 10 ; Average over these many swap steps
threshold = 1 ; Do not swap if < threshold
cyl0-r = 0.8 ;
cyl0-up = 0.5 ; Split cylinder 0 upper extension (nm)
cyl0-down = 0.5 ; Split cylinder 0 lower extension (nm)
cyl1-r = 0.8 ; same for other channel
cyl1-up = 0.5
cyl1-down = 0.5
I’m not sure exactly what the problem is, but here are a few things I’d check or try.
Firstly, what is the charge difference between your compartments, could it be too big?
Second, is your double membrane system stable without using the CompEl protocol?
Thirdly, are your compartments well defined? The CompEl output file should give you the z position of both compartment boundaries, are they always well inside the membrane?
Thanks for answering. I have used gmx genions to distribute the ions at the beginning, hence I don’t think the charge difference would be that great… also with visual inspection, I can see that the ion is roughly uniform distributed.
The system run when it single layer membrane and double layer membrane, and was able to run 1 microsecond without crashing… I am think whether this is a problem of force field too… because this system is set up with martini3 forcefield…
I have set up the CompEI compatment boundaries as the COM of the protein channel, as it was suggested in the tutorial. The protein still within the membrane, i.e. inbetween the lipid system with visual inspection and confirmed with density plot… Therefore, I really have no idea what is wrong…
at the very start of your simulation, i.e. at the first time step, are there swaps done? Could you paste the start of your swapions.xvg output file? Maybe that sheds some light on it.
There is no swapions.xvg output file, however this is the new error:
SWAP: Warning! Step 150, ion 167909 moved from Domain_B to Domain_A
but did not pass cyl0 or cyl1 as defined in the .mdp file.
Do you have an ion somewhere within the membrane?
Currently because I am randomly insert the ions, I am hoping this is the error of the system isn’t equalibrate yet… I am going to do a 1 microsecond simulation and see whether this error go away.
By visual inspection, it seems like the ions are going through the channel and not anywhere else.
there must be some .xvg file from the computational e-phys module, maybe you renamed it by using the -deffnm option of mdrun? What command line did you use for calling mdrun?
It could be that you have a small pocket of water somewhere within the membrane (could be just a couple of water molecules). If so, I would suggest to remove these manually before any more simulations, as in my experience it can take ages for the water to find a way out of the membrane. When filling in ions randomly, it could happen that one of the waters within the membrane is replaced by an ion, causing the trouble. Ions in the membrane will cause exactly the Warning that you see - they will move between both compartments without crossing a channel.
Whether and how many water/ion exchanges the CompEl protocol executes depends on the number of ions in both compartments. With ions in the membrane, these numbers cannot be computed reliably, so it’s best to avoid such a scenario.
When I look at the .xvg output file, I notice the following things:
Could it be that there is an enormous difference in charge (Delta q) between the compartments? It looks like the simulation starts with Delta q = 78z and then continues at Delta q = 2042z after the first swap. This will probably lead to a voltage that is more than the membrane can handle. But maybe I am wrong, the numbers are hard to find as the formatting is messed up. There should be a line in the .xvg output file starting with # Requested charge imbalance is ... that you have not pasted here - could you look that up please?
The number of water molecules in the two compartments is very different (80k vs 50k) - is this intentional? Normally I would set up the system by placing two identical membrane systems on top of each other, which would result in the same number of water molecules in both compartments at the beginning of the simulation (and also would yield a Delta q = 0z to start with).
The .xvg output file tells you the index numbers of the problematic ions (e.g. ion 167909) Can you check in Pymol where this ion is? This might shed some light on the problem.
Thanks for your input and knowledge, so currently I am going to remove all the water from the membrane and you are correct that the current charge different is: # Requested charge imbalance is Q(A)-Q(B) = -2042 e. We have such a big charge difference is based on experimental results, i.e. the K intracellular is 110 mM and Extracellular is 4 mM.
In the experiment, they provided an external voltage onto the membrane to open and close position of the protein. However in simulation, doing so required adding external electric force onto the membrane. The current consenses is we are not going to provide an electric field in the system, where we are going just let the K ion flow naturally. (We are not sure whether this will be a good approach either, but we welcome some input here.)
We put the membrane where one is mirror of the other protein, hence it will form a outer and inner membrane region. I originally put is in the 50 nm box where the 2membrane is at 15 and 40 nm, respectively. No idea what happened there during initiation of the simulation…
But 100 mM versus 4 mM K only refers to the concentrations of K+ ions. For the resulting charge difference you need to take into account also the other ions like Cl-.
Yes but in your case a Q(A)-Q(B) = -2042 e will result in a strong transmembrane electric field. If you don’t want that, you need to set up a simulation with Q(A)-Q(B) = 0 e. You could estimate the transmembrane potential using the gmx potential tool.
After using your suggestion and removing all the water molecules and ion atoms from from membrane, the error still occurs:
Do you have an ion somewhere within the membrane?
SWAP: Warning! Step 577800, ion 164136 moved from Domain_B to Domain_A
but did not pass cyl0 or cyl1 as defined in the .mdp file.
Do you have an ion somewhere within the membrane?
However, using the VMD to visualise the 164136 ion, it is actually inside the channel, but diffusing slowly (it is doing in and out motion)… henceforth, is there another way to make this warning to go away?
I have also double the cylinder size in the mdp file:
the split cylinder should be large enough in radius so that an ion has to move through it to get to the other membrane side (it does not hurt if the radius is somewhat larger than the pore). On top of that, if you test for ion/water exchanges every N steps, then the height of the cylinder should be larger than what an ion can move in N steps (so that no ion can slip through unnoticed by the protocol.) To test for ions slipping through, you can set N to 10 or even 1, just to test whether the warning goes away with that.
Thanks for the quick response. I have set the coupl-steps = 1, however, the same problem still occurs. The height of the channel is much greater than what an ion can move. From VMD, I can see that the ions is in the channel but not moving out. I am assuming due to that the ion not transporting into the inner membrane (or from Zone A to Zone B). I think this may be the cause of that error? Therefore, anything else that I can do?
So this is in my .mdp file right now:
swapcoords = Z ; Swap positions: no, X, Y, Z
swap-frequency = 100 ; Swap attempt frequency
split-group0 = layer1 ; Defines compartment boundary
split-group1 = layer2 ; Defines other compartment boundary
massw-split0 = yes ; use mass-weighted center?
massw-split1 = yes
solvent-group = W ; Group containing the solvent molecules
iontypes = 2 ; Number of different ion types to control
iontype0-name = K
iontype0-in-A = 48 ; concentraion of 0.114 mol/l, inner:outer membrane 110:4 K ion ratio
iontype0-in-B = 1322
iontype1-name = CL
iontype1-in-A = -1
iontype1-in-B = -1
coupl-steps = 1 ; Average over these many swap steps
threshold = 1 ; Do not swap if < threshold
cyl0-r = 1.6 ; Split cylinder 0 radius (nm);
cyl0-up = 1.5 ; Split cylinder 0 upper extension (nm)
cyl0-down = 1.5 ; Split cylinder 0 lower extension (nm)
cyl1-r = 1.6 ; same for other channel
cyl1-up = 1.5
cyl1-down = 1.5
Thank you for your advice. I fixed the problem by exclude the protrude protein from the membrane in the file. Now I am trying to understand the swapions.xvg file, as such I don’t understand what does the s1, s2… etc. mean:
Which I am trying to understand, and I believe it is incorrect due to it has an expontential decrease of potential… Therefore, is this due to I fixed the number of Chlorine ion in the zone A and zone B, which caused this negative potential?
the s0, s1, etc are labels for the data columns. You can use xmgrace for example to quickly display the content of an .xvg file graphically.
For gmx potential, be sure to a) include the whole system, as all charges contribute to the z-averaged potential (if you just use a part of all atoms, the output will be meaningless), b) pass the -correct flag to gmx potential indicating that the potential on both upper and lower end of the periodic box can be assumed to be identical (therefore the system must have zero total charge).