//This function was implemented by Sudar Purushothanan
double ComputeGC(void)
{
      MOLECULE *tempNode;
      MOLECULE *tempDegreeNode;
      MOLECULE *neighborNode;
      MOLECULE *thirdNode;
      double gc=0;
      bool prim=false;
      double tempQ;
      double tempZ;
      int primaryQuadsFound=0;
      int secondaryQuadsFound=0;
      int secondaryDegrees=0;
      int primaryDegrees=0;
 
      tempNode=MoleculesList;
      while (tempNode)
      {
            for (int i=0;i< tempNode->linksCount;i++)
            {
                  tempQ=0;
                  tempZ=0;
                  neighborNode=GetNodeBasedOnNumber(tempNode->linksTo[i]);
                  if (neighborNode)
                  {
                        for (int j=0;j<neighborNode->linksCount;j++)
                        {     
                              thirdNode=GetNodeBasedOnNumber(neighborNode->linksTo[j]);
                              if ((!thirdNode) || (thirdNode->number==tempNode->number))
                              {
                              }
                              else
                              {
                                    for (int k=0;k<thirdNode->linksCount;k++)
                                    {
                                          for (int l=0;l<tempNode->linksCount;l++)
                                          {
                                                if (thirdNode->linksTo[k]==neighborNode->number)
                                                {
                                                }
                                                else
                                                {
                                                      if (thirdNode->linksTo[k]==tempNode->linksTo[l])
                                                      {
                                                            for (int m=0;m<thirdNode->linksCount;m++)
                                                            {
                                                                  if (thirdNode->linksTo[m]==tempNode->number)
                                                                  {
                                                                        primaryQuadsFound++;
                                                                        prim=true;
                                                                  }
                                                            }
                                                            if (!prim)
                                                            {
                                                                  secondaryQuadsFound++;
                                                            }
                                                            prim=false;
                                                      }
                                                }
                                          }
                                    }
                              }
                        }
                  }
            }
            primaryDegrees=tempNode->linksCount;
            for (int n=0;n<tempNode->linksCount;n++)
            {
                  tempDegreeNode=GetNodeBasedOnNumber(tempNode->linksTo[n]);
                  if (tempDegreeNode)
                  {
                        secondaryDegrees+=tempDegreeNode->linksCount-1;
                  }
                  else
                  {
                        secondaryDegrees+= 0;
                  }
            }
            tempQ=(double)(primaryQuadsFound+secondaryQuadsFound);
            tempQ/=2;
            tempZ=(double)(primaryDegrees*(primaryDegrees-1)*(primaryDegrees-2))/2;
            tempZ+=(double)(secondaryDegrees*primaryDegrees*(primaryDegrees-1))/2;
            if (tempQ>tempZ)
            {
                  tempQ=tempQ;
                  tempZ=tempZ;
            }
            if (tempZ>0)
            {
                  gc+=(double)(tempQ/tempZ);
                  tempNode->gc = (double)(tempQ/tempZ);
            }
            else
            {
                  tempNode->gc = 0;
            }
            tempNode=tempNode->next;
            primaryQuadsFound=0;
            secondaryQuadsFound=0;
            primaryDegrees=0;
            secondaryDegrees=0;
      }
      double kval=double(NumberOfInteractions);
      kval/=double(NumberOfMolecules);
      gc=(double)(gc/NumberOfMolecules);
      return gc;
}