Here is some trips and tricks based on our calculations with Linux, Gromacs and Plumed. You may find it useful in your research (but it is just made for my personal convenience)

1. How to generate stripped PDB with deleted ANISOU

grep -v ANISOU my-pdb-with-ANISOU.pdb > my-pdb-without-ANISOU.pdb

2. Position restraint procedure for C-alphas

gmx make_ndx -f confout.part0001.gro -o index.ndx

gmx genrestr -f confout.part0001.gro -n index.ndx -fc 500 500 500 -o posreca.itp

3. Trajectory handling

Processing trajectory

#If you have jump problem during visualisation

gmx trjcat -f traj_comp.part*.xtc -o trajall.xtc

echo 1 | gmx trjconv -s ~/pfs/HIV/ref.tpr -f trajall.xtc -pbc nojump -o trajnoj.xtc

echo 3 1 | gmx trjconv -s ~/pfs/HIV/ref.tpr -f trajnoj.xtc -fit rot+trans -o trajfit.xtc

gmx trjcat -f trajfit.xtc -dt 10 -o tfit.xtc

#If there is problem with periodic boundary conditions in Gromacs smulation

trjcat -f traj*.part*.xtc -o trajall.xtc

echo 1 | trjconv -s ../gpwtest.tpr -f trajall.xtc -pbc mol -o trajmol.xtc

echo 3 1 | trjconv -s ../gpwtest.tpr -f trajmol.xtc -fit rot+trans -o trajfit.xtc

trjcat -f trajfit.xtc -dt 10 -o tfit.xtc

Extract a single frame for restart

gmx trjconv -s plm.tpr -f traj.part0001.trr -dump 1800 -o conf1800.gro

Concatenate trajectories

gmx trjcat -f traj.part*.xtc -o trajall.xtc

Combining trajectories from different simulations

trjcat -f R1/traj_r1_times.xtc R2/traj_r2_times.xtc -cat -nosort -o combine.xtc

the flags:

-cat : Do not discard double time frames

-nosort:  Sort trajectory files (not frames)

To make the output trajectory (combine.xtc) more organized we can use the -timestep flag

trjcat -f R1/traj_r1_times.xtc R2/traj_r2_times.xtc -cat -nosort -o combine.xtc -timestep 1

-timestep :  Change time step between input frames (ps)

Fix PBC problem (try one of them):

gmx trjconv -s ref.tpr -f trajall.xtc -pbc mol -o trajmol.xtc

gmx trjconv -s ref.tpr -f trajall.xtc -pbc nojump -o trajnoj.xtc

Align trajectory to the reference frame

trjconv -s ref.tpr -f trajmol.xtc -fit rot+trans -o trajfit.xtc

Calculating distance

Use "gmx distance -n index.ndx -f traj.xtc -s topol.tpr -oav -oall" with index file write in this way:

[CA_1 – CA_2]

300 402

Sparse out trajectory (write every 10 ps) or cut out certain interval

trjcat -f trajfit.xtc -dt 10 -o tfit.xtc

trjcat -f trajfit.xtc -b 100 -e 140 -o trajfit_interestingpart.xtc

Fix the volume for NVT simulation (after NPT equilibration)

gmx eneconv -f ener.part*.edr -o enerall.edr

gmx energy -f enerall.edr -o vol.xvg

xmgrace vol.xvg

head vol.xvg

grep '^....0' vol.xvg

grep '^...00' vol.xvg

(just a trick to see only points represented in the trr trajectory, every 10 or 100 ps)

Select a frame (near the end) that has a volume close to the average of the last 10 nanoseconds or so. Dump it as above.

Covariance analysis, projections on principal components (modes), and visualization of modes

gmx covar -s plm.tpr -f trajfit.xtc

gmx anaeig -s plm.tpr -v eigenvec.trr -f trajfit.xtc -proj proj.xvg -last 10

gmx anaeig -s plm.tpr -v eigenvec.trr -f trajfit.xtc -first 4 -last 4 -extr ext4.pdb -nframes 100

(gives an artificial trajectory between the extreme points of mode 4, use Rock in VMD to see it rocking)

#convert trajectory to xyz format, examples:

writexyz_delfirst ../protstart.pdb tfit.xtc

writexyz gpwtest.tpr ../tfit.xtc

writexyz -parm7 ../com_solvated.top -crdbox ../test.mdcrd

In case of .nc or netcdf file (.nc or netcdf format is binary whereas .mdcrd is ascii). In that case

writexyz -parm7 ../com_solvated.top -nc ../md_dry.nc

(uses VMD, so takes the same arguments as vmd, outfile is called out.xyz)

Convert AMBER trajectory to GROMACS or vice-versa or other general trajectory conversion

install mdtraj

$ sudo pip install mdtraj

if pip, numpy and cython is not pre-installed you need to install these packages.
sudo apt-get install python-pip
sudo pip install numpy cython
then execute
sudo pip install mdtraj

Upon installation you can exceute the conversion using the following command or similar arguments.
/usr/local/bin/mdconvert ~/Erik/plasmepsin_correct/R2/md_total.nc -o traj_r2.xtc

Extend Simulation with Plumed

#gmx convert-tpr extends the simulations by the amount of time in ps which is mentioned with the -extend flag. -s = input .tpr file -o = output .tpr file

gmx convert-tpr -s gpwtest.tpr -o gpwtest.tpr -extend 30000

Clustering analysis with G_cluster

g_cluster is capable clustering using RMSD.

/home/s/sbhakat/pfs/bin/gromacsbin462/bin/g_cluster_plumed_mpi_d -f ../../Test_combine/combine_all_times.xtc -s ../../prot0.gro -n ../../../1D_opt/index.ndx -xvg none -wcl 5 -cutoff 0.04 -dt 100 -method gromos -minstruct 150 -clid clust-id.xvg

4. Calculate free-energy surface after a simulation

One-dimensional:

plumed sum_hills --hills 'HILLS' --outfile fes --stride 100 --bin 100

bin = number of bins

stride = how often it outputs intermediate files (time between files = stride*pace*timestep)

Two-dimensional:

plumed sum_hills --hills HILLS --outfile fes --stride 500 --bin 100,100

Two-dimensional FES projected on one CV:

plumed sum_hills --hills HILLS --outfile fes --stride 500 --bin 100,100 --idw d1 --kt 2.494

Note : The result depends on the temperature (2.494 = kT with T=300K)

5. Plotting examples

(use “head COLVAR” to see what the columns mean)

gpu 1:2 COLVAR

gpu '1:($2*10)' COLVAR   (note the quotes to protect $ from bash)

gpuwp 1:3 COLVAR

gpsev 3 COLVAR    (assumes one x column, plots the next 3 columns)

gp fes*.dat     (plots several files in the same graph)

Several separate plots in the same figure:

gnuplot> set multiplot layout 4,2

multiplot> sp 'fes2.dat' u 1:2:3

multiplot> sp 'fes3.dat' u 1:2:3

multiplot> sp 'fes4.dat' u 1:2:3

multiplot> sp 'fes5.dat' u 1:2:3

multiplot> sp 'fes6.dat' u 1:2:3

multiplot> sp 'fes7.dat' u 1:2:3

multiplot> sp 'fes8.dat' u 1:2:3

multiplot> sp 'fes9.dat' u 1:2:3

multiplot> unset multiplot

Plot two-dimensional FES:

gnuplot> set pm3d map

gnuplot> sp ‘fes1.dat’ u 1:2:3

Histogram

#Usage: histocol.py file column [nbin [ start stop] ]

histocol.py COLVAR 2 100 > hist1.dat

Running averages

#Usage: runav.py file windowlength

runav.py COLVAR 100 > COLVAR_av100

(averages all columns except the x column)

Average and standard deviation

nohead 1 COLVAR | head -500 | awk '{print $2}' |stdev

Test switching function (and plot histogram in the same plot)

p 'hist6.dat' u ($1*0.1):(0.003*$2) w l, (1-((x-d0)/r0)**n)/(1-((x-d0)/r0)**m), d0=0.29,n=6,m=12,r0=0.05

Install Matlab from an ISO image

mkdir Mathworks

sudo mount -o loop -t iso9660 ~/matlab/R2015a-glnxa64.iso Mathworks/ (give system password if asked)

cd Mathworks

./install

umount Mathworks

Unmount a directory

You can use the following command

fusermount -u /some/directory/

Here are the options of fusermount

fusermount: [options] mountpoint
Options:
 -h print help
 -V print version
 -o opt[,opt...] mount options
 -u unmount
 -q quiet
 -z lazy unmount

A more detailed source to read regarding umount is here

Downloading big files from Cluster scp stalls problem

rsync -avz --progress sbhakat@aurora.lunarc.lu.se:/lunarc/nobackup/users/sbhakat/BACE/Metadynamics/3TPL/COM_TIP3P_DIST/trajfit.xtc Aurora_download/BACE/Metadynamics/3TPL/

File conversion:

moving the files starting with # to a new file

mv "#traj.part0001.xtc.1#" traj.part0000.xtc


Automatic job submission script for Gromacs

#!/bin/sh
#
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#
# Job name
#SBATCH -J HIV_gromacs
#
# No. of nodes and no. of processors per node
#SBATCH -N 1
#SBATCH --tasks-per-node=16
#SBATCH --exclusive
#
# Time needed to complete the job
#SBATCH -t 48:00:00
#
module add gcc/4.8.1
module add openmpi/1.8.1/gcc/4.8.1
module add cuda/6.5
if [ -e state.cpt  ]; then
  extra='-cpi state.cpt'
fi
# GMX-5.0.4
#
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#
mpirun -npernode 16 /home/per/lun/gromacs5bin/bin/mdrun_mpi -s hiv300ns.tpr $extra -ntomp  1 -gpu_id 0000000011111111 -v -maxh 48
rc=$?
echo rc $rc
if [ $rc != 0 ]; then echo >stop; fi


if [ -f stop  ]; then
  exit
fi
if [ -f confout.*  ]; then
  exit
fi

sleep 10m
sbatch job.sh

To compare tICA output from MSMBuilder (RawpositionFeaturiser) with Plumed Driver use the TYPE=SIMPLE as TYPE=OPTIMAL although is perfectly linear it is doing some kind of normalization inside Plumed.

PCAVARS REFERENCE=reference_corrected.pdb TYPE=SIMPLE LABEL=tica
PCAVARS REFERENCE=reference_corrected.pdb TYPE=OPTIMAL LABEL=ticaopt
PRINT ARG=tica.* FILE=COLVAR
PRINT ARG=ticaopt.* FILE=COLVAROPT

Task: Run a job over VPN even if VPN disconnects

nohup command   &

Task: Numbering the DESRES trajectory

awk -v CONVFMT=%.5f '{a=a+0.00025; $1=a; print $0}' COLVAR_hbondnohead >COLVAR

Task: Find the max of a particular column

sort -nrk1,1 filename | head -1 | cut -d ' ' -f3 #this is for getting max for column 1

Replacing ‘blank space’ by ‘comma’

tr ' ' ',' <tmp.txt >output.txt

A rule of thumb bandwidth estimator for histograms https://en.wikipedia.org/wiki/Kernel_density_estimation

Count number of columns

awk '{print NF}' COLVAR_sincos | sort -nu | tail -n 1