From S.R.Kilvington@soton.ac.uk  Mon Jul  4 10:25:47 1994
Received: from mail.soton.ac.uk  for S.R.Kilvington@soton.ac.uk
	by www.ccl.net (8.6.9/930601.1506) id JAA04700; Mon, 4 Jul 1994 09:42:51 -0400
Received: from localhost (srk@localhost) by mail.soton.ac.uk (8.6.4/2.12) id OAA01790 for CHEMISTRY@ccl.net; Mon, 4 Jul 1994 14:42:03 +0100
From: Simon Kilvington <S.R.Kilvington@soton.ac.uk>
Message-Id: <199407041342.OAA01790@mail.soton.ac.uk>
Date: Mon, 4 Jul 94 14:42:02 BST
To: CHEMISTRY@ccl.net
Subject: molecular superpositioning program
X-Mailer: ELM [version 2.3 PL11]


Dear net people,

        With the recent interest in molecular super-positioning I
thought I'd post a copy of my program that I sent to Mike Smith (who
originally wrote to the list wanting a program to overlay two
molecules).

	The program takes two PDB files and two lists of atom numbers
and produces a translated/rotated version of the second PDB file so the
specified atoms are overlayed as closely as possible.

        I'm sending the source with this message, hopefully it should
get through uncorrupted. If there are any problems due to more than 80
characters per line (my text editor is set up with about 150 columns)
I could ftp it to you or something.

        There are 3 source files, "overlay.c", "utils.c" and "utils.h".
After you've chopped them out of this file, make the program by doing

cc -o overlay overlay.c utils.c

	You'll also need a "-lm" flag on SGI's as for some reason the
maths part of the standard C library isn't automatically linked in.

        Then to overlay the two PDB files you do

overlay <pdb1> <atom_list1> <pdb2> <atom_list2>

        This will produce a file called "<pdb2>.super". And tell you
the RMS distance between the atoms you specified to be overlayed.

        The atom list files should have one atom id per line. The
atom id can be just its index number, or you can specify the chain
too. The format of each line is...

<atom index>[<white space><chain id>]

        where things in brackets [] are optional.

eg.
118
4560
etc

or

118     A
118     B
4560    B
etc


        Obviously there must be the same number of atoms in each list.

        I hope this is useful. If there are any problems or queries send me a
message.

        yours,

	Simon Kilvington (srk@uk.ac.soton)

The source...

=================================================================
=================================================================

/*
   overlay.c

   Simon Kilvington, 1994
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "utils.h"

#define MAXFILENAMELEN		256
#define PDBLINELEN		82		/* length of a line in a pdb file, inc space for "\n\0" terminator */
#define INDEXCOLUMN		6		/* position of atom number in ATOM records */
#define CHAINCOLUMN		21		/* position of chain id in ATOM records */
#define XPOSCOLUMN		30		/* position of x coord in ATOM records */
#define YPOSCOLUMN		38
#define ZPOSCOLUMN		46
#define COORDLEN		8		/* width of coord fields */
#define COORDDPS		3		/* number of decimal places for coords */
#define SUPERPOSTOLERANCE	0.000000175	/* when superpositioning give up when tan(rotation needed) is less than this, ~tan(0.00001degs) */
#define MAXSUPERPOSITS		100000		/* or we need more than this number of iterations */

typedef struct
{
   char   pdb[PDBLINELEN];
   int    index;				/* pdb atom number */
   char   chain;				/* chain id, resolves atom number clashes */
   vector pos;
} atomdetails;

typedef struct
{
   int         natoms;
   atomdetails atom[1];				/* array extends to atom[natoms-1] */
} fragdetails;

#define FRAGHDRSIZE	(sizeof(fragdetails) - sizeof(atomdetails))

fragdetails *readpdbdata(char *);
void writepdbdata(char *, fragdetails *);

int *readatomlist(char *, fragdetails *, int *);
int findatomindex(fragdetails *, int, char);

void superpositionfrags(int, fragdetails *, int *, fragdetails *, int *);
void calcsuperposmatrix(int, fragdetails *, int *, fragdetails *, int *, matrix *, vector *, vector *);

/*
   main
   params should be <pdb file1> <atom file1> <pdb file2> <atom file2>
   the <atom files> should be lists of atoms (and, optionally, chain id's) to be overlapped, there must be the same number in each file
   produces a pdb file called "<pdb file2>.super"
*/

int
main(int argc, char **argv)
{
   fragdetails *pdb1, *pdb2;
   int *atom1, *atom2, natoms1, natoms2, a;
   float rms;
   vector dr;
   char filename[MAXFILENAMELEN];

   if(argc != 5)
   {
      printf("Useage: overlay <pdb1> <atom list1> <pdb2> <atom list2>\n");
      return 1;
   }

/* read in the 	files */
   if((pdb1 = readpdbdata(argv[1])) == NULL)
      return 1;
   if((atom1 = readatomlist(argv[2], pdb1, &natoms1)) == NULL)
      return 1;
   if((pdb2 = readpdbdata(argv[3])) == NULL)
      return 1;
   if((atom2 = readatomlist(argv[4], pdb2, &natoms2)) == NULL)
      return 1;
   if(natoms1 != natoms2)
   {
      printf("Both atom lists must have the same number of atoms\n");
      return 1;
   }

   printf("Moving \"%s\"\n", argv[3]);
   superpositionfrags(natoms1, pdb1, atom1, pdb2, atom2);

   strcpy(filename, argv[3]);
   strcat(filename, ".super");
   writepdbdata(filename, pdb2);

/* report the rms distance between atom1[] and atom2[] */
   rms = 0.0;
   for(a=0; a<natoms1; a++)
      rms += vector_len2(vector_sub(&pdb1->atom[atom1[a]].pos, &pdb2->atom[atom2[a]].pos, &dr));
   rms = sqrt(rms / natoms1);
   printf("RMS distance between atoms to overlap: %.3f\n", rms);
   printf("Number of atom pairs: %d\n", natoms1);

   free(pdb1);
   free(atom1);
   free(pdb2);
   free(atom2);

   return 0;
}

fragdetails *
readpdbdata(char *filename)
{
   fragdetails *frag, *newptr;
   atomdetails *aptr;
   FILE *file;
   char buffer[PDBLINELEN];
   BOOL err;

   if((file = fopen(filename, "r")) == NULL)
   {
      printf("Unable to read \"%s\"\n", filename);
      return NULL;
   }

   printf("Reading \"%s\"\n", filename);

   frag = malloc(FRAGHDRSIZE);
   frag->natoms = 0;

   err = FALSE;
   while(!feof(file) && !err)
   {
   /* read in the line */
      fgets(buffer, PDBLINELEN, file);
   /* if this is an ATOM record add it to the growing fragment */
      if(strncmp(buffer, "ATOM  ", 6) != 0)
	 continue;
      frag->natoms++;
      newptr = realloc(frag, FRAGHDRSIZE + (frag->natoms * sizeof(atomdetails)));
      if(err = (newptr == NULL))
	 continue;
      frag = newptr;
      aptr = &frag->atom[frag->natoms-1];
      strncpy(aptr->pdb, buffer, PDBLINELEN);
      aptr->index = atoi(&buffer[INDEXCOLUMN]);
      aptr->chain = buffer[CHAINCOLUMN];
      aptr->pos.x = atof(&buffer[XPOSCOLUMN]);
      aptr->pos.y = atof(&buffer[YPOSCOLUMN]);
      aptr->pos.z = atof(&buffer[ZPOSCOLUMN]);
   }

   fclose(file);

   if(err)
   {
      printf("No memory to read \"%s\"\n", filename);
      free(frag);
   }

   return (err) ? NULL : frag;
}

void
writepdbdata(char *filename, fragdetails *frag)
{
   FILE *file;
   int a;
   atomdetails *aptr;

   if((file = fopen(filename, "w")) == NULL)
   {
      printf("Unable to write to \"%s\"\n", filename);
      return;
   }

   printf("Writing \"%s\"\n", filename);

   for(a=0; a<frag->natoms; a++)
   {
      aptr = &frag->atom[a];
      sprintf(&aptr->pdb[XPOSCOLUMN], "%*.*f", COORDLEN, COORDDPS, aptr->pos.x);
      sprintf(&aptr->pdb[YPOSCOLUMN], "%*.*f", COORDLEN, COORDDPS, aptr->pos.y);
      sprintf(&aptr->pdb[ZPOSCOLUMN], "%*.*f", COORDLEN, COORDDPS, aptr->pos.z);
      fprintf(file, "%s\n", aptr->pdb);
   }

   fclose(file);

   return;
}

int *
readatomlist(char *filename, fragdetails *frag, int *natoms)
{
   int *list, *newptr, index, bpos;
   char chain, buffer[PDBLINELEN];
   FILE *file;
   BOOL err;

   *natoms = 0;

   if((file = fopen(filename, "r")) == NULL)
   {
      printf("Unable to read \"%s\"\n", filename);
      return NULL;
   }

   printf("Reading \"%s\"\n", filename);

   list = malloc(sizeof(int));	/* so we can use realloc from the start */

   err = FALSE;
   while(!feof(file) && !err)
   {
      (*natoms)++;
      newptr = realloc(list, (*natoms) * sizeof(int));
      if(err = (newptr == NULL))
      {
	 printf("No memory read \"%s\"\n", filename);
	 continue;
      }
      list = newptr;
      fgets(buffer, PDBLINELEN, file);
      sscanf(buffer, "%d%n", &index, &bpos);
      while(buffer[bpos] != '\0' && isspace(buffer[bpos]))
	  bpos++;
      chain = (buffer[bpos] == '\0') ? ' ' : buffer[bpos];
      if(err = ((list[(*natoms)-1] = findatomindex(frag, index, chain)) == -1))
	 printf("Error while reading \"%s\"; PDB atom with index %d%c does not exist\n", filename, index, chain);
   }

   fclose(file);

   if(err)
      free(list);

   return (err) ? NULL : list;
}

int
findatomindex(fragdetails *frag, int index, char chain)
{
   int atom;
   BOOL match;

   atom = 0;
   match = FALSE;
   while(!match && atom < frag->natoms)
   {
      match = (frag->atom[atom].index == index && frag->atom[atom].chain == chain);
      if(!match)
         atom++;
   }

   return (atom == frag->natoms) ? -1 : atom;
}

/*
   superpositionfrags
   this moves yfrag so the specified atoms in it and xfrag overlap as closely as possible
   the xatom and yatom arrays should be natoms long and hold the numbers of the atoms that are required to be overlappped
   uses the Ferro-Hermans algorithm, Acta Cryst. 33 (1977) p345-347
*/

void
superpositionfrags(int natoms, fragdetails *xfrag, int *xatom, fragdetails *yfrag, int *yatom)
{
   matrix spos, ypos, result;
   vector cgx, cgy;
   int a;
   atomdetails *aptr;

   matrix_alloc(&spos, 3, 3, TRUE);
   matrix_alloc(&ypos, 3, 1, TRUE);

   calcsuperposmatrix(natoms, xfrag, xatom, yfrag, yatom, &spos, &cgx, &cgy);

   for(a=0; a<yfrag->natoms; a++)
   {
      aptr = &yfrag->atom[a];
      matrix_setelement(&ypos, 0, 0, aptr->pos.x - cgy.x);
      matrix_setelement(&ypos, 1, 0, aptr->pos.y - cgy.y);
      matrix_setelement(&ypos, 2, 0, aptr->pos.z - cgy.z);
      matrix_mult(&spos, &ypos, &result);
      aptr->pos.x = cgx.x + matrix_getelement(&result, 0, 0);
      aptr->pos.y = cgx.y + matrix_getelement(&result, 1, 0);
      aptr->pos.z = cgx.z + matrix_getelement(&result, 2, 0);
      matrix_forget(&result);
   }

   matrix_forget(&ypos);
   matrix_forget(&spos);

   return;
}

/*
   calcsuperposmatrix
   as superpositionfrags, but yfrag is unmoved by the experience
   on entry spos should be a matrix_alloc'd 3x3 matrix
   on exit spos is the superposition matrix, cgx is the centre of gravity of xatom, and cgy is the centre of gravity of yatom
   to superposition the two frags yfrag should be moved to "cgx + spos(y - cgy)" where y are the original coords of yfrag
   !!! this routine depends on the fact that a "vector" is composed of 3 consecutive floats that represent x, y and z coords !!!
*/

#define veccoord(VEC, AXIS)	(*(((float *) (VEC)) + AXIS))
#define xcoord(ATOM, AXIS)	veccoord(&(xfrag->atom[xatom[ATOM]].pos), AXIS)
#define ycoord(ATOM, AXIS)	veccoord(&(yfrag->atom[yatom[ATOM]].pos), AXIS)

enum {XAXIS, YAXIS, ZAXIS};

void
calcsuperposmatrix(int natoms, fragdetails *xfrag, int *xatom, fragdetails *yfrag, int *yatom, matrix *spos, vector *cgx, vector *cgy)
{
   matrix corr;
   int i, j, k, p, q, r, itno;
   float oneovern, sum, sigma, gamma, dist, qk, rk;
   BOOL rotated;

   if(natoms < 3)
   {
      printf("Not enough atoms (%d) to superposition molecules\n", natoms);
      return;
   }

   matrix_alloc(&corr, 3, 3, TRUE);

/* calc cgx and cgy */
   cgx->x = cgx->y = cgx->z = 0.0;
   cgy->x = cgy->y = cgy->z = 0.0;
   for(i=0; i<natoms; i++)
   {
      vector_add(cgx, &xfrag->atom[xatom[i]].pos, cgx);
      vector_add(cgy, &yfrag->atom[yatom[i]].pos, cgy);
   }
   oneovern = 1.0 / (float) natoms;
   vector_scale(cgx, oneovern);
   vector_scale(cgy, oneovern);

/* initialise the correlation matrix, c'' in the paper */
   for(j=XAXIS; j<=ZAXIS; j++)
   {
      for(k=XAXIS; k<=ZAXIS; k++)
      {
	 sum = 0.0;
	 for(i=0; i<natoms; i++)
	    sum += (xcoord(i, j) - veccoord(cgx, j)) * (ycoord(i, k) - veccoord(cgy, k));
	 matrix_setelement(&corr, k, j, sum);
      }
   }
/* spos is a bit easier to do, this is M in the paper */
   matrix_makeunit(spos);

   itno = 0;
   do
   {
   /* calc the rotation needed about each axis */
      rotated = FALSE;
      for(p=XAXIS; p<=ZAXIS; p++)
      {
	 q = (p+1 <= ZAXIS) ? p+1 : XAXIS;
	 r = (q+1 <= ZAXIS) ? q+1 : XAXIS;
      /* calc the rotation needed, angle = arctan(sigma/gamma) */
	 sigma = matrix_getelement(&corr, r, q) - matrix_getelement(&corr, q, r);
	 gamma = matrix_getelement(&corr, q, q) + matrix_getelement(&corr, r, r);
      /* if the angle is not too small update the two matrices */
	 if(fabs(sigma / gamma) > SUPERPOSTOLERANCE)
	 {
	    dist = (float) sqrt(gamma*gamma + sigma*sigma);
	    if(dist != 0.0)
	    {
	       for(k=XAXIS; k<=ZAXIS; k++)
	       {
	       /* the superposition matrix */
		  qk = matrix_getelement(spos, q, k);
		  rk = matrix_getelement(spos, r, k);
		  matrix_setelement(spos, q, k, (gamma*qk + sigma*rk) / dist);
		  matrix_setelement(spos, r, k, (-sigma*qk + gamma*rk) / dist);
	       /* and the correlation matrix */
		  qk = matrix_getelement(&corr, q, k);
		  rk = matrix_getelement(&corr, r, k);
		  matrix_setelement(&corr, q, k, (gamma*qk + sigma*rk) / dist);
		  matrix_setelement(&corr, r, k, (-sigma*qk + gamma*rk) / dist);
	       }
	    }
	    rotated = TRUE;			/* make a note of the fact that a rotation was still needed */
	 }
      }
      itno++;
   }
   while(rotated && itno < MAXSUPERPOSITS);	/* until we didnt rotate about any axis, or we give up */

   matrix_forget(&corr);

   return;
}

#undef veccoord
#undef xcoord
#undef ycoord

=================================================================
=================================================================

/*
   utils.c

   Simon Kilvington, 1994
*/

#include "utils.h"

float
vector_len2(vector *vp)
{
   return (vp->x * vp->x) + (vp->y * vp->y) + (vp->z * vp->z);
}

vector *
vector_add(vector *vp1, vector *vp2, vector *vp3)
{
   vp3->x = vp1->x + vp2->x;
   vp3->y = vp1->y + vp2->y;
   vp3->z = vp1->z + vp2->z;

   return vp3;
}

vector *
vector_sub(vector *vp1, vector *vp2, vector *vp3)
{
   vp3->x =vp1->x - vp2->x;
   vp3->y =vp1->y - vp2->y;
   vp3->z =vp1->z - vp2->z;

   return vp3;
}

vector *
vector_scale(vector *v, float scalar)
{
   v->x *= scalar;
   v->y *= scalar;
   v->z *= scalar;

   return v;
}

BOOL
matrix_alloc(matrix *m, int rows, int cols, BOOL quit)
{
   if(!(m->element = (matrix_element_t *) malloc(rows * cols * sizeof(matrix_element_t))))
   {
      printf("Out of memory in matrix_alloc\n");
      if(quit)
	 exit(1);
      return FALSE;
   }

   m->rows = rows;
   m->cols = cols;

   return TRUE;
}

matrix *
matrix_makeunit(matrix *m)
{
   int i, j;

   for(i=0; i<m->rows; i++)
      for(j=0; j<m->cols; j++)
	 matrix_setelement(m, i, j, (i==j) ? 1.0 : 0.0);

   return m;
}

matrix *
matrix_mult(matrix *m1, matrix *m2, matrix *r)
{
   int i, j, k;
   matrix_element_t e;

   if(m1->cols != m2->rows)
   {
      printf("Matrices are of incompatible sizes in matrix_mult\n");
      exit(1);
   }

   matrix_alloc(r, m1->rows, m2->cols, TRUE);

   for(i=0; i<m1->rows; i++)
   {
      for(j=0; j<m2->cols; j++)
      {
         e=0;
         for(k=0; k<m1->cols; k++)
            e += matrix_getelement(m1, i, k) * matrix_getelement(m2, k, j);
         matrix_setelement(r, i, j, e);
      }
   }

   return r;
}

=================================================================
=================================================================

/*
   utils.h

   Simon Kilvington, 1994
*/

#include <math.h>
#include <stdio.h>

typedef enum {FALSE, TRUE} BOOL;

typedef struct
{
   float x, y, z;
} vector;

float vector_len2(vector *);
#define vector_len(VP)		sqrt(vector_len2(VP))

vector *vector_add(vector *, vector *, vector *);
vector *vector_sub(vector *, vector *, vector *);

vector *vector_scale(vector *, float);

typedef float matrix_element_t;

/* a matrix is just a dope vector and a ptr to an array of elements */
typedef struct
{
   int rows, cols;
   matrix_element_t *element;
} matrix;

#define matrix_elementptr(mat, i, j)		(((mat)->element) + ((j) * (mat)->rows) + (i))
#define matrix_setelement(mat, i, j, val)	(*(matrix_elementptr(mat, i, j)) = (val))
#define matrix_getelement(mat, i, j)		(*(matrix_elementptr(mat, i, j)))

BOOL matrix_alloc(matrix *, int, int, BOOL);
#define matrix_forget(mat)			(free((mat)->element), (mat)->rows = (mat)->cols = 0)

matrix *matrix_makeunit(matrix *);

matrix *matrix_mult(matrix *, matrix *, matrix *);

=================================================================
=================================================================

Here endeth todays source.