Analysing crossing waters in membrane with LOOS

[T98G]

Hi,

I'm trying to use the crossing-waters tool to analyse the number of water molecules that cross a membrane in a trajectory from gromacs. I've centralized my membrane with gmx trjconv, and devided my trajectory into 100ns chunks

crossing-waters npt.gro com_centered_0_100.xtc 10 20 'resname == "POPG"'

But no matter what I do I get the same absurd result:

# Total frames = 501
# Number of waters = 13056
#AtomID	Lifetime	Entered	Exited	ExitedPositive

I've centered the membrane in all possible ways, but nothing seems to work, this result is exactly the same for different trajectory chunks, like the last or the first 100 ns of a 2µs simulation.

My system has 17920 water molecules and it's just impossible for 13056 water molecules to cross the membrane in 100 ns. In all publications regarding water pores the cumulative flux of water is in the order of 2500 molecules over a period of 500 ns. So it's just impossible this is right. The result's change to 256 if I use 'name == "P"' for the selection, but they don't change whether I use a 100 ns chunk of trajectory or the whole 2µs trajectory.

I really don't know what I'm doing wrong.

Thanks

[agrossfield]

Hi,

I'll do my best to help you.

That output indicates that it thinks there are 13056 waters in the system and none of them crossed during the simulation. However, from your description, that's clearly not what you expected (the zero crossings in 500 frames could be right, but the number of waters isn't).

I'm reasonably certain the problem is that the selection string at the end is supposed to be the water, and you're giving it one of the lipid species. What happens if you leave the selection string out entirely? There's a default water selection built into the code, which often works but isn't universal. If it doesn't give you there expected 17920, you might need to customize the selection -- we guessed about common namings for water, but can't match everything. If it doesn't work, you might try something like

'resname == "WAT" && !hydrogen'

replacing "WAT" with whatever your water residue name is. I'm assuming you're using a standard 3-site water model, and you'd need something a little different for a 4-site model.

Separately, if you're going to use LOOS more than once it's worth running gmxdump2pdb.pl to get a fake psf to use instead of a gro file:

gmx dump -s foo.tpr | gmxdump2pdb.ps PREFIX --constraints

where you'd replace foo.tpr with your system's tpr file, and PREFIX with whatever you'd like the resulting pdb and psf files to be called.

LOOS is willing to use gro files, but they have almost no metadata in them -- this conversion to psf gives you a system file with bonds, charges, masses, residue names, and a guess at segment names, which makes a lot of analysis tasks simpler. As an added benefit, VMD can do more with the PSF than it can with the gro (e.g. if it's a coarse-grained simulation, the psf will let it draw bonds correctly).

Just to be clear: using the gro isn't what's causing this problem. It's just that going forward, you're better off using a richer file format.

Last, you might want to take a look at some of the user tutorials for LOOS, available at https://github.com/GrossfieldLab/loos/wiki/Tutorials-for-Users

In particular, there's a tutorial focused on getting started with LOOS for gromacs simulations.

I hope this helps -- please write back if it doesn't, or you have more questions.

Cheers,

Alan

[T98G]

HI, thank you so much for you reply, it was incredibly helpful

I had missunderstood that the selection was suposed to be the membrane, I just tried 'resname == "SOL" && name == "OW"' for my selection and it seems to have worked.

I got an output like this (with many more lines):

70172 4 200 204 1
70202 1 203 204 1
70325 1 203 204 1

But I don't quite understant this data, the number of water molecules that crossed is in the second or in the last collumn?

Thank you very much,

Tarcillo

[agrossfield]

Hi Tarcillo,

A good starting place for almost any LOOS tool is to run with the fullhelp option: crossing-waters --fullhelp
This gives a detailed description of the expected input, output, usually something about the algorithm, references if any, and potential issues to be aware of.

Along the same lines, the first couple of lines of output of almost every tool contains the command-line used to invoke the program and some kind of documentation of the output. For example:

# crossing-waters 'system.psf' 'traj.dcd' '10' '20' - alan (Tue Mar 13 14:32:49 2012) {/directory/you/were/working/in} [2.0.0 120313]
# Total frames = 719
#AtomID	Lifetime	Entered	Exited	ExitedPositive
38546	2	0	2	-1
25136	2	4	6	1

would be the header and 1st two lines of data. The lines marked by "#" are comments, showing you the header -- the last one gives the meaning of the columns of data.

Each data line is a single water that crossed. The first number is identifies which water it is (we use the atom ID from the file). The second is the amount of time the water spent inside the membrane, in frames (e.g. if your trajectory has frames 1 ns apart, these would be 2 ns). The next 2 columns are the frame number where the molecule entered and exited the membrane, and the last number is whether it exited at negative or positive z.

By itself, this output isn't super-helpful -- you'll need to post-process it to get what you want. If you just want the overall rate of crossing, just count the number of lines and divide by the length of the trajectory. If you want to know whether the rate is changing, you can group by when the water entered (or left) and compute mini-rates (e.g. how many waters crossed during the first 10 ns, the second 10 ns, etc). Or, you might want to know how long waters tend to stay in the membrane, in which case you'd focus on the lifetime column. There are enough possibilities that I didn't want to hard-code them in (not to mention, when I wrote it you'd be lucky to run long enough to get 10 or 20 events total, so it didn't seem worth it).

Hope this helps,

Alan

[T98G]

Thank you very much for your reply, it was very helpful. I'm working with a membrane containing a water pore, so there are lots of water crossings, but I found that loos crossing-waters is perfect for me, as I didn't want to spend the time to write such tool myself.

Thanks a lot,

best regards,

Tarcillo