Commit 424c6b22 authored by lheinz's avatar lheinz
Browse files

changed gitignore

parent a7318505
Installation of g_permute
-------------------------
*********************
PREREQUISITE: GROMACS 2018
*********************
g_permute consists of two parts, the 'lap' library in 'g_permute/liblap' and the tool
itself, 'g_permute' in 'g_permute/src.
Both parts rely on the GROMACS MD package for some functions.
The code has been tested with and requires GROMACS 2018
Adjust Makefile to your GROMACS installation (compiled with shared library support)
------------------------------------------------------
1) cd into 'src' and edit the Makefile
2) set PREFIX to the directory you want the final library to go
Make sure the path you enter here is in your LD_LIBRARY_PATH when running g_permute
3) specifiy GMXDIR, GMXLIB and GMXINC to point to your Gromacs installation
make sure you have compiled gromacs with shared libraries enabled
(i.e. with 'configure --enable-shared')
4) in src/, type 'make', then 'make install'
Have fun !
User Manual Version 1.1
#######################
1. What is g_permute ?
======================
g_permute is a tool to relabel the solvent molecules in a molecular
dynamics trajectory such that their distance to a reference position
becomes minimal. The result is a new trajectory, where the solvent
molecules are centered around their reference positions, rather than
exploiting the full configuration space as in the original trajectory.
This renders established entropy estimation methods applicable and makes
diffusive systems accessible to improved fit functions. For further details
see the publications cited on the webpage mentioned below.
g_permute was written in the computational biophysics group (Helmut Grubmuller)
at the MPI for biophysical chemistry in Goettingen (Germany)
and is available under the download section
http://www.mpibpc.mpg.de/home/grubmueller/
Feel free to contact juergenhaas@gmx.net for further information.
2. Installation
===============
The prerequisite GROMACS package and its dependencies can be installed
automatically with a supplied script.
Please consult the ./INSTALL file for details!
3. Usage
========
3.1 Input arguments
-------------------
g_permute takes as input arguments
three filenames
-f the input trajectory
-n an index file
-r the reference structure
and one numerical argument
-m the size of the solvent molecules
its output goes into a file specified by
-o the output trajectory
-f can be in any Gromacs-readable format
-r specifies the reference structure. In the course of g_permute, the
solvent molecules will be relabeled in each frame of the input trajectory
such that the distance to the solvent of this reference structure
becomes minimal.
The reference structure thus has to contain the same topology as any
frame of the input trajectory. Again, any Gromacs-readable format is
possible.
If none given the positions in the first frame are used as reference frame.
-n is an index file. It has to contain
1) one index group specifying all solvent atoms. The atom indices have to
be grouped into molecules. e.g, for TIP4P water (four atoms, one oxygen,
two hydrogens and a dummy atom), the first molecule corresponds to
indices 1-4, the second to 5-8, ...
2) one index group specifying one key atom in each solvent molecule.
The distance to the reference structure is computed as
the total distance of these key atoms to their positions
in the reference structure
The user is responsible to make sure that each solvent
molecule of group 1 contains only one key atom of group 2
--- example for TIP4P water ---
[ oxygen_atoms ]
1 5 9 13 etc.
[ solvent_atoms ]
1 2 3 4 5 6 etc.
In the case of relabeling TIP4P water, the first
group contains the indices of all O atoms of the water molecules. The
second group contains the indices of all atoms.
-m The size of a molecule (e.g., 4 for TIP4P water). The solvent group
is internally split into blocks of size 'm', which explains the
grouping condition mentioned above.
After startup, g_permute will ask you to choose from your index file
- the group for the distance calculation
(i.e., group 2) from above, for TIP4P water: the oxygens)
- the solvent group (group 1 from above)
3.2 optional switches
---------------------
g_permute takes up to 5 optional switches, where the switches -b and -e are to be used for frames selection, while -com and -rm_pbc affect the distance calculation itself. The switch -dt reduces the number of frames to include only those where time t modulo the timestep of the trajectory dt is equal the time indicated.
-com center of mass (COM) mode instead of default picking of first atom (oxygen in case of water as solvent)
-rm_pbc calculate distance using periodic boundary conditions, default is off
-b Begin at the specified frame in ps
-e End at the specified frame in ps
-dt Only consider frames when time t MOD dt == first time (ps), otherwise disregard frame.
3.3 Data output
---------------
Output is written to a file specified by -o, there are two optional output files containing COM coordinates, either permuted or in original order.
-o the output trajectory
-oref trajectory containing COM coordinates in original order (optional)
-orp trajectory containing the permuted COM coordinates (optional)
4. Example run
==============
The simplest command to permute TIP4P water molecules (4 atoms per molecule) in a
GROMACS based trajectory is:
g permute -m 4 -s start_struc.pdb -f traj.xtc -o permute.xtc -n index.ndx
g permute then reads in the topologies of the system from start struc.pdb
and the trajectory from traj.xtc, referring to the groups listed on index.ndx. The
first group will be used for the determining the distances and the second one for the solvent group.
The permuted trajectory will be written to permute.xtc. Note that the output
trajectory can also be written directly to a multi-model pdb file.
5. Hints and practicalities
===========================
5.1 Composing index files
-------------------------
Since we cannot do without an index file, we recommend to use the GROMACS
tool make_ndx. A call to make_ndx requires a structure file in 'pdb' or
a run input file in the GROMACS 'tpr' format.
make_ndx -f start_struc.pdb
On a box filled exclusively with tip4p water molecules make_ndx will give
the following options:
There are: 216 OTHER residues
There are: 0 PROTEIN residues
There are: 0 DNA residues
Analysing Other...
0 System : 864 atoms
1 SOL : 864 atoms
nr : group ! 'name' nr name 'splitch' nr Enter: list groups
'a': atom & 'del' nr 'splitres' nr 'l': list residues
't': atom type | 'keep' nr 'splitat' nr 'h': help
'r': residue 'res' nr 'chain' char
"name": group 'case': case sensitive 'q': save and quit
Choose now all oxygen atoms:
>a OW
Found 216 atoms with name OW
2 OW : 216 atoms
press the <ENTER_KEY> to list the groups
0 System : 864 atoms
1 SOL : 864 atoms
2 OW : 216 atoms
nr : group ! 'name' nr name 'splitch' nr Enter: list groups
'a': atom & 'del' nr 'splitres' nr 'l': list residues
't': atom type | 'keep' nr 'splitat' nr 'h': help
'r': residue 'res' nr 'chain' char
"name": group 'case': case sensitive 'q': save and quit
Choose to save to 'index.ndx' and quit.
>q
You now should have the newly created index.ndx.
5.2 Visualisation of trajectories to investigate effect of permutation
----------------------------------------------------------------------
>>>visualise with vmd:
vmd start_struc.pdb permute.xtc
select graphical representations
Enter into Selection field 'resname OW and resid 1 to 10'
Draw Style > drawing method > VdW
Draw Style > colouring > resid
Trajectory > Draw multiple frames > 0:100
>>>repeat this for the unpermuted trajectory:
vmd start_struc.pdb traj.xtc
In the window 'VMD Main' click on the 'D' to de/active the display of that
individual trajectory to compare the two trajectories. In the permuted trajectory
the oxygen atoms will be centered around a reference position and thus cover a
smaller area than in the original trajectory.
5.3 Memory issues
-----------------
The number of solvent molecules contained in these groups is limited by
the available system memory (RAM), so if you encounter messages like not
enough memory, play with the number of molecules in the solvent group
(SOL) in the case of water.
5.4 dash / sh compatibility issues
----------------------------------
On ubuntu systems you might encounter shell errors, which can be overcome by explicitely running 'bash g_permute/instgmx331.sh' or 'g_permute/instgmx331.sh'.
/************************************************************************
*
* gnrl.h
version 1.0 - 21 june 1996
author Roy Jonker, MagicLogic Optimization Inc.
header file for general stuff
*
**************************************************************************/
/*************** CONSTANTS *******************/
#if !defined TRUE
#define TRUE 1
#endif
#if !defined FALSE
#define FALSE 0
#endif
/*************** DATA TYPES *******************/
typedef int boolean;
/*
* This file contains some functions to simplify memory handling
* of matrices and vectors. It was taken from
* the Jonker / Volgenant LAP code and adapted slightly
* see http://www.magiclogic.com/assignment.html for details
*/
#include <memory.h>
#include "lap.h"
#include <stdio.h>
#include <stdlib.h>
void terminate(void);
int *alloc_ints(int len);
void free_ints(int *ints);
cost *alloc_costs(int len);
void free_costs(cost *costs);
cost **alloc_costmatrix(int width, int height);
void free_costmatrix(cost **cmat, int width, int height);
/************************************************************************
*
* lap.h
version 1.0 - 21 june 1996
author Roy Jonker, MagicLogic Optimization Inc.
header file for LAP
*
**************************************************************************/
//#include <typedefs.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#ifndef __LAP_H
#define __LAP_H
#ifdef CPP
#define EXPORT extern "C"
#endif
#ifndef CPP
#define EXPORT extern
#endif
/*************** CONSTANTS *******************/
#define BIG 100000
/*************** TYPES *******************/
typedef int row;
typedef int col;
typedef double cost;
/*************** FUNCTIONS *******************/
EXPORT cost lap(int dim, cost **assigncost, int *rowsol, int *colsol, cost *u, cost *v);
EXPORT void checklap(int dim, cost **assigncost, int *rowsol, int *colsol, cost *u, cost *v);
#endif /* __LAP_H */
Licensing information for liblap
--------------------------------
This library is an implementation of an algorithm proposed in
R. Jonker and A. Volgenant (University of Amsterdam)
"A Shortest Augmenting Path Algorithm for Dense and Sparse Linear Assignment Problems", Computing 38, 325-340 (1987)
The implementation distributed with g_permute is a slightly changed
version of the original implementation of the authors which can be found at
http://www.magiclogic.com/assignment.html
Note that the code is *not* distributed under GPL. Instead, it is
licensed as follows:
"The code is Copyright 2003 MagicLogic Systems Inc., Canada.
These codes are made available in the hope that they will be useful,
but without any warranty. Also, we [the authors] cannot accept any
responsibility for their use.
The codes are available only for non-commercial use. Please contact
MagicLogic for commercial application."
With permission of the authors, this code has been slightly changed
for use with g_permute.
/************************************************************************
*
* lap.cpp
version 1.0 - 4 September 1996
author: Roy Jonker @ MagicLogic Optimization Inc.
e-mail: roy_jonker@magiclogic.com
Code for Linear Assignment Problem, according to
"A Shortest Augmenting Path Algorithm for Dense and Sparse Linear
Assignment Problems," Computing 38, 325-340, 1987
by
R. Jonker and A. Volgenant, University of Amsterdam.
*
*************************************************************************/
//#include "smalloc.h"
#include "gromacs/utility/smalloc.h"
#include "gnrl.h"
#include "lap.h"
#include <math.h>
#define VERYSMALL 0.00001
double lap(int dim,
double **assigncost,
col *rowsol,
row *colsol,
double *u,
double *v)
/* input: */
/* dim - problem size */
/* assigncost - cost matrix */
/* output: */
/* rowsol - column assigned to row in solution */
/* colsol - row assigned to column in solution */
/* u - dual variables, row reduction numbers */
/* v - dual variables, column reduction numbers */
{
boolean unassignedfound;
row i, imin, numfree = 0, prvnumfree, f, i0, k, freerow;
col j, j1, j2 = 0, endofpath = 0, last = 0, low, up;
double min = 0, h, umin, usubmin, v2, lapcost;
int loopcnt;
row *unassigned = NULL; /* list of unassigned rows. */
col *collist = NULL; /* list of columns to be scanned in various ways. */
col *matches = NULL; /* counts how many times a row could be assigned. */
double *d = NULL; /* 'cost-distance' in augmenting path calculation. */
row *pred = NULL; /* row-predecessor of column in augmenting/alternating path. */
snew(unassigned, dim);
snew(collist, dim);
snew(matches, dim);
snew(d, dim);
snew(pred, dim);
#ifdef _VERBOSE
for(i=0; i<dim; i++){
for(j=0; j<dim; j++)
printf("%f ", assigncost[i][j]);
printf("\n");
}
i = j = 0;
#endif
/* init how many times a row will be assigned in the column reduction. */
for (i = 0; i < dim; i++)
matches[i] = 0;
/* COLUMN REDUCTION */
for (j = dim-1; j >= 0; j--) /* reverse order gives better results. */
{
/* find minimum cost over rows. */
min = assigncost[0][j];
imin = 0;
for (i = 1; i < dim; i++)
if (assigncost[i][j] < min)
{
min = assigncost[i][j];
imin = i;
}
v[j] = min;
if (++matches[imin] == 1)
{
/* init assignment if minimum row assigned for first time. */
rowsol[imin] = j;
colsol[j] = imin;
}
else
colsol[j] = -1; /* row already assigned, column not assigned. */
}
/* REDUCTION TRANSFER */
for (i = 0; i < dim; i++)
if (matches[i] == 0) /* fill list of unassigned 'free' rows. */
unassigned[numfree++] = i;
else
if (matches[i] == 1)/* transfer reduction from rows that are assigned once. */
{
j1 = rowsol[i];
min = BIG;
for (j = 0; j < dim; j++)
if (j != j1)
if (assigncost[i][j] - v[j] < min)
min = assigncost[i][j] - v[j];
v[j1] = v[j1] - min;
}
/* AUGMENTING ROW REDUCTION */
loopcnt = 0; /* do-loop to be done twice. */
do
{
loopcnt++;
/* scan all free rows. */
/* in some cases, a free row may be replaced with another one to be scanned next. */
k = 0;
prvnumfree = numfree;
numfree = 0;/* start list of rows still free after augmenting row reduction. */
while (k < prvnumfree)
{
i = unassigned[k];
k++;
/* find minimum and second minimum reduced cost over columns. */
umin = assigncost[i][0] - v[0];
j1 = 0;
usubmin = BIG;
for (j = 1; j < dim; j++)
{
h = assigncost[i][j] - v[j];
if (h < usubmin) {
if (h >= umin)
{
usubmin = h;
j2 = j;
}
else
{
usubmin = umin;
umin = h;
j2 = j1;
j1 = j;
}
}
}
i0 = colsol[j1];
if (umin < usubmin) /* - VERYSMALL) */
/* change the reduction of the minimum column to increase the minimum */
/* reduced cost in the row to the subminimum. */
v[j1] = v[j1] - (usubmin - umin);
else /* minimum and subminimum equal. */
if (i0 >= 0) /* minimum column j1 is assigned. */
{
/* swap columns j1 and j2, as j2 may be unassigned. */
j1 = j2;
i0 = colsol[j2];
}
/* (re-)assign i to j1, possibly de-assigning an i0. */
rowsol[i] = j1;
colsol[j1] = i;
if (i0 >= 0) {/* minimum column j1 assigned earlier. */
if (umin < usubmin)/* - VERYSMALL) */
/* put in current k, and go back to that k. */
/* continue augmenting path i - j1 with i0. */
unassigned[--k] = i0;
else
/* no further augmenting reduction possible. */
/* store i0 in list of free rows for next phase. */
unassigned[numfree++] = i0; }
}
}
while (loopcnt < 2); /* repeat once. */
/* AUGMENT SOLUTION for each free row. */
for (f = 0; f < numfree; f++)
{
freerow = unassigned[f];/* start row of augmenting path. */
/* Dijkstra shortest path algorithm. */
/* runs until unassigned column added to shortest path tree. */
for (j = 0; j < dim; j++)
{
d[j] = assigncost[freerow][j] - v[j];
pred[j] = freerow;
collist[j] = j; /* init column list. */
}
low = 0; /* columns in 0..low-1 are ready, now none. */
up = 0; /* columns in low..up-1 are to be scanned for current minimum, now none. */
/* columns in up..dim-1 are to be considered later to find new minimum, */
/* at this stage the list simply contains all columns */
unassignedfound = FALSE;
do
{
if (up == low) /* no more columns to be scanned for current minimum. */
{
last = low - 1;
/* scan columns for up..dim-1 to find all indices for which new minimum occurs. */
/* store these indices between low..up-1 (increasing up). */
min = d[collist[up++]];
for (k = up; k < dim; k++)
{
j = collist[k];
h = d[j];
if (h <= min)
{
if (h < min) /* new minimum. */
{
up = low; /* restart list at index low. */
min = h;
}
/* new index with same minimum, put on undex up, and extend list. */
collist[k] = collist[up];
collist[up++] = j;
}
}
/* check if any of the minimum columns happens to be unassigned. */
/* if so, we have an augmenting path right away. */
for (k = low; k < up; k++)
if (colsol[collist[k]] < 0)
{
endofpath = collist[k];
unassignedfound = TRUE;
break;
}
}
if (!unassignedfound)
{
/* update 'distances' between freerow and all unscanned columns, via next scanned column. */
j1 = collist[low];
low++;
i = colsol[j1];
h = assigncost[i][j1] - v[j1] - min;
for (k = up; k < dim; k++)
{
j = collist[k];
v2 = assigncost[i][j] - v[j] - h;
if (v2 < d[j])
{
pred[j] = i;
if (v2 == min) {/* new column found at same minimum value */
if (colsol[j] < 0)
{
/* if unassigned, shortest augmenting path is complete. */
endofpath = j;
unassignedfound = TRUE;
break;
}
/* else add to list to be scanned right away. */
else
{
collist[k] = collist[up];
collist[up++] = j;
}}
d[j] = v2;
}
}
}
}
while (!unassignedfound);