[gmx-users] Membrane Simulation in Gromacs - Extracting waters from hydrophobic region

[Antoniel A. S. Gomes]

Hi Gromacs users.

Our lab is working on membrane simulations using gromacs. The main problem to build a system is extract waters from aliphatic carbons. Gromacs website have a script called keepbyz.sh (link: (link: http://www.gromacs.org/Documentation/How-tos/Membrane_Simulations) but it have some problems.

If your system have 10000 atoms or more, the second and third columns become one and the script changes the axis search from Z to Y. This happens because the script gets the line's value from left to right. Based on this, we created a new one to avoid that problem.

Benefits?

- It dispenses all the text pre- and post-treatment (just needs the solvate .gro file as input) and generates only one output file, ready to be opened in vmd program;
- It can be used in .gro files with or without velocities;
- It shows how many waters are present in the new file (useful to edit the .top file and proceed the simulation);
- It is two times faster.

Our script is set to select the Z axis. If you want to get the Y axis, just change the $(NF) value to $(NF-1) and $(NF-4) in lines 12 and 23, respectively. Same idea to X axis (changing to 2 and 5). The editable variables for upper and lower threshold are shown as up and down (lines 2 and 3), respectively. The script is presented below:

############

#!/bin/bash
up=5.414
down=1.788
watercount=0

read_n() { for i in $(seq 3); do read || return; echo "$REPLY"; done; }

cat $1 | grep -v 'SOL' | sed \$d >solvate_clean.gro
cat $1 | grep 'SOL' | (while lines=$(read_n); do
vel=$(echo "$lines" | head -1 | awk '{print NF}')
	if [ $vel -le "6" ]; then
	atom=$(echo "$lines" | grep ' OW' | awk '{print $(NF)}')
	exp1=$(echo "$atom >= $up" | bc)
	exp2=$(echo "$atom <= $down" | bc)
		if [ $exp1 == 1 ] || [ $exp2 == 1 ] ;then
			echo "Water Kept"
			echo "$lines">>solvate_clean.gro
			watercount=$(expr $watercount + 1)
		else
			echo "Water Discarded"
		fi
	elif [ $vel -gt "6" ]; then
	atom=$(echo "$lines" | grep ' OW' | awk '{print $(NF-3)}')
	exp1=$(echo "$atom >= $up" | bc)
	exp2=$(echo "$atom <= $down" | bc)
		if [ $exp1 == 1 ] || [ $exp2 == 1 ] ;then
			echo "Water Kept"
			echo "$lines">>solvate_clean.gro
                        watercount=$(expr $watercount + 1)
		else
			echo "Water Discarded"
		fi
	fi
done
echo -e "\nYour new system have $watercount waters. Edit your .top file\n")
tail -1 $1 >>solvate_clean.gro
sed -i '2s/.*/'$(expr $(cat solvate_clean.gro | wc -l) - 3)'/' solvate_clean.gro

############

Usage: ./cleanwater.sh file.gro

Feel free to use and share it if you want.

Best,

Antoniel A. S. Gomes

Bacharel em Ciências Biológicas, UEPB - Campus V
Mestre em Biologia Celular e Molecular, UFPB - Campus I
Doutorando em Biologia Geral e Aplicada, UNESP - Instituto de Biociências - Botucatu
Lattes: http://lattes.cnpq.br/2330777925652141