/******************************************************************/
/* Calculate Rg (radius of gyration) from ATOM coordinates of PDB */
/*      M. Kojima, Tokyo University of Pharmacy and Life Sciences */
/*                                                2018.01/25      */
/* 1. Save this file as calcRg.c                                  */
/* 2. Compile and build executable file such as calcRg.exe        */
/*  $ gcc -o calcRg.exe calcRg.c                                  */
/* 3. Run the program                                             */
/*  $ calcRg.exe PDBfilename                                      */  
/******************************************************************/   
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define LBMAX 81          /* size of line buffer */

typedef struct
{
        char name[5];
        float coord[3];
	float mass;
} atom_data;

atom_data *ATOM;
int num_atoms;
char record[5];

void quit()
{
	free(ATOM); 
        ATOM = NULL;
}
	
void ReadPDB(char *file_name)
{
        FILE *fp;
	char line[LBMAX];
	int i = 0;

	if((fp = fopen(file_name, "r")) == NULL) {
	  fprintf(stderr, "Could not open %s\n", file_name);
  	  exit(1);
        }
	else
	  while(fgets(line, LBMAX, fp) != NULL) {
	   sscanf(line,"%4c",record);
	   if(strcmp(record,"ATOM") == 0){

	     if(i == 0) { 	
	       ATOM = (atom_data *)malloc(sizeof(atom_data));
	     } else {	  
	       ATOM = (atom_data *)realloc(ATOM, (i+1)*sizeof(atom_data));
	     } 

	     sscanf(line, "%*12c%4c%*14c%8f%8f%8f", ATOM[i].name, &ATOM[i].coord[0],  &ATOM[i].coord[1], &ATOM[i].coord[2]);

	     switch(ATOM[i].name[1]) {
             case 'H':
               ATOM[i].mass = 1.0; 
	       break;
	     case 'C':
	       ATOM[i].mass = 12.0;
	       break;
	     case 'N':
	       ATOM[i].mass = 14.0;
	       break;
	     case 'O':
      	       ATOM[i].mass = 16.0;
	       break;
	     case 'S': 
      	       ATOM[i].mass = 32.0;
	       break;
	     default:
	       ATOM[i].mass = 0.0;
               printf("Illegal atom name: %4s\n", ATOM[i].name);
	       break;
	     }	  
	     i++;
	  } 
        }
        num_atoms = i;
        printf("Total number of atoms: %d\n",num_atoms);	
        fclose(fp);
}	

float rg()
{
	float center[3] = {0.0, 0.0, 0.0};
	float total_mass = 0.0;
        float rg = 0.0;
	int i, j;

	for(i = 0; i < num_atoms; i++)
	{
	  for(j = 0; j < 3; j++)
	  {
	    center[j] += ATOM[i].mass*ATOM[i].coord[j];
	  }
	  total_mass += ATOM[i].mass;
	}

	for(j = 0; j < 3; j++)
	{
	  center[j] /= total_mass;
  	}

	printf("Total mass: %6.2f Da\n", total_mass);

	for(i = 0; i < num_atoms; i++)
	{
	  for(j = 0; j < 3; j++)
	  {
	    rg += (ATOM[i].coord[j] - center[j])*(ATOM[i].coord[j] - center[j])*ATOM[i].mass;
	  }
	}

        return sqrt(rg / total_mass);
}

int main(int argc, char *argv[])
{
	int i;
	
	if (argc < 2) {
	  fprintf(stderr, "usage: calcRg.exe PDBfilename \n");
	  exit(1);
	}
	
	ReadPDB(argv[1]);
	printf("Rg= %8.3f A\n", rg());
        quit();

	return 0;
}