/*
************************************************************************
*
*   ss.C -
*
*   Copyright (c) 1995
*
*   ETH Zuerich
*   Institut fuer Molekularbiologie und Biophysik
*   ETH-Hoenggerberg
*   CH-8093 Zuerich
*
*   All Rights Reserved
*
*   Date of last modification : 95/03/07
*   Pathname of SCCS file     : /export/home3/cb/garant-5.0/src/SCCS/s.ss.C
*   SCCS identification       : 1.2
*
************************************************************************
*/
/**************************************************************************/
/*                            ss.cc                                       */
/*                                                                        */
/* class to load and generate spin systems                                */
/**************************************************************************/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "bib.h"		/* library with bonds and atoms              */
#include "log.h"		/* general log file                          */
#include "map.h"		/* classes to compare graphs                 */
#include "ss.h"

/* Structure *****************************************************************/

int Structure::fragAdd(char *fn, int type)
{
   BibAccess bib;
   int res;

   res = bib.fragLoad(fn,type);
   res &= fragLoad();
   res &= atomLoad();
   res &= bondLoad();

   return res;
}

int Structure::atomLoad()
{
   AtomRequest Ar;
   Atom *A,*A1,*A2;
   Iterator it1,it2;
   Frag *fr;
   int res = 1;
   PseudRequest pR;
   LiAtom liPseud;
   int oldNr=-1;

   Ar.start();

   while(!Ar.eof) {
      fr = (Frag *)Frags.find(Frag(Ar.nrFrag));
      A = new Atom(Ar.nr,Ar.name,Ar.type,Wert(Ar.shift,Ar.delshift),Ar.prob,
                   fr);
      if(!Atoms.insert(A)) {
         delete A;
      }
      Ar.read();
   }

   pR.start();
   while(!pR.eof) {
      if(oldNr != pR.pseudnr) {
         oldNr = pR.pseudnr;
         A1=(Atom *)liPseud.start(it1);
         while(A1!=0) {
            A2=(Atom *)liPseud.start(it2);
            while(A2!=0) {
               A1->pseud.insert(A2);
               A2=(Atom *)liPseud.next(it2);
            }
            A1=(Atom *)liPseud.next(it1);
         }
         liPseud.remove();
      }
      A=(Atom *)Atoms.find(Atom(pR.atnr));
      if(A) liPseud.insert(A);
      pR.read();
   }

   return(res);
}

int Structure::fragLoad()
{
   FragRequest Fr;
   Frag *F;
   int res = 1;
   static int isfirstDest=1;

   Fr.start();

   while(!Fr.eof) {
      F = new Frag(Fr.type, Fr.nrInt, Fr.nrExt, Fr.name, Fr.equivNr);
      if(!Frags.insert(F)) {
         delete F;
      } else {
         if(Fr.type==DESTFRAG) {
            if(isfirstDest) {
               isfirstDest=0;
               Frags.firstDest=Fr.nrInt;
               Frags.lastDest=Fr.nrInt;
            } else {
               if(Frags.firstDest>Fr.nrInt) Frags.firstDest=Fr.nrInt;
               if(Frags.lastDest<Fr.nrInt) Frags.lastDest=Fr.nrInt;
            }
         }
      }
      Fr.read();
   }
   return(res);
}

int Structure::bondLoad()
{
   BondRequest Br;
   int res = 1;

   Br.start();
   while(!Br.eof) {
      res = bond(Br.nr1, Br.nr2);
      if(res) res = bond(Br.nr2, Br.nr1);
      Br.read();
   }
   if(res) return 1;
   else return 0;
}

int Structure::bond(int a1,int a2)
{
   Atom *A1,*A2;
   Bond *B;
   int res = 1;

   A1= (Atom *)Atoms.find(Atom(a1));
   A2= (Atom *)Atoms.find(Atom(a2));

   B = new Bond(A1,A2);

   if(!A1->rels.insert(B)) {
      res = 0;
      delete B;
   }
   if(res == 0) return(0);
   else return(a2);
}

void Structure::listAt(char *t)
{
   Log log;

   log.w("\n... list of Atoms\n");
   Atoms.list(6,t);
}

void Structure::loadCoord(char *fn, int fragtype)
{
   CoordRequest cR;
   int cnt = 0;
   char str[MAXLINE];
   Log log;
   Atom *at;
   Coord *coord;
   OrdLiAtom ordAts=OrdLiAtom(Atoms);
   static int idStruct=0;

   cR.start(fn);
   while(!cR.eof) {
      at = ordAts.find(cR.nrExt,cR.name,fragtype);
      if(at!= 0) {
         coord = new Coord(cR.coord,idStruct);
         at->coords.insert(coord);
	 cnt++;
      } else {
         sprintf(str,"\n... WARNING: unknown atom %s %d",cR.name,cR.nrExt);
         log.warn(-1,str);
      }
      cR.read();
   }
   idStruct++;
   sprintf(str,"\n... %d coordinates defined",cnt);
   log.w(str);
}

void Ss::jMake()
{
   JTopol T;
   int res;
   int number=0;
   Log log;
   char str[MAXLINE];

   res=T.start();
   while(res) {
      number += T.match(this);
      T.cleanup();
      res=T.next();
   }

   sprintf(str,
       "\n... %d scalar couplings generated",number);
   log.w(str);
}

void Ss::statNOEMake(int f,int l)
{
   StatNOETopol T(f,l);
   int res;
   int number=0;
   Log log;
   char str[MAXLINE];

   res=T.start();
   while(res) {
      number += T.match(this);
      T.cleanup();
      res=T.next();
   }
   sprintf(str,
       "\n... %d dipolar couplings generated",number);
   log.w(str);
}

void Ss::distNOEMake(int limitOk, float limit, float p, int usefloat)
{
   int number=0;
   Log log;
   char str[MAXLINE];
   Atom *a1,*a2,*a1p,*a2p;
   Coord *c1,*c2,*cdummy;
   Iterator i1,i2,i1p,i2p,idummy;
   float dist,diff,minDist=0.0;
   int i;
   int nrOk,ok;
   Noe *noe;

   limit = limit * limit;
   a1 = (Atom *)Atoms.start(i1);
   while(a1 != 0) {
      i2 = i1;
      a2 = (Atom *)Atoms.next(i2);
      while(a2 != 0 && a1->proton()) {
        nrOk=0;
        cdummy=(Coord *)a1->coords.start(idummy);
        while(cdummy!=0 && a2->proton() && a1!=a2) {
           ok=0; minDist=-1.0;
           a1p=(Atom *)a1->pseud.start(i1p);
           if(a1p==0 || !usefloat)a1p=a1;
           while(a1p!=0) {
              a2p=(Atom *)a2->pseud.start(i2p);
              if(a2p==0 || !usefloat)a2p=a2;
              while(a2p!=0) {
                 if(a1p != a2p &&
                    a1p->frag->type == a2p->frag->type &&
                    a1p->p != 1.0 && a2p->p != 1.0) {
                    c1=(Coord *)a1p->coords.get(idummy);
                    c2=(Coord *)a2p->coords.get(idummy);
                    if(c1!=0 && c2!=0) {
                       dist = 0.0;
                       for(i=0;i<3;i++) {
                          diff = fabs(c1->coord[i] - c2->coord[i]);
                          diff -= (a1p->pseudoCorr + a2p->pseudoCorr);
                          if(diff<0.0) diff=0.0;
                          dist += diff * diff;
                       }
                       if(dist <= limit) ok++;
                       if(minDist<0.0 || minDist>sqrt(dist)) minDist=sqrt(dist);
                    }
                 }
                 if(a2p->pseud.getn()>0 && usefloat)
                    a2p=(Atom *)a2->pseud.next(i2p);
                 else a2p=0;
              }
              if(a1p->pseud.getn()>0 && usefloat)
                 a1p=(Atom *)a1->pseud.next(i1p);
              else a1p=0;
           }
           if(ok>0)nrOk++;
           cdummy=(Coord *)a1->coords.next(idummy);
        }
        if(nrOk >= limitOk) {
           if(!a1->rels.Index::find(Noe(a1,a2))) {
              noe = new Noe(a1,a2);
              noe->prob = p;
              noe->minDist = minDist;
              (void)a1->rels.insert(noe);
              number++;
           }
           if(!a1->rels.Index::find(Noe(a2,a1))) {
              noe = new Noe(a2,a1);
              noe->prob = p;
              noe->minDist = minDist;
              (void)a2->rels.insert(noe);
              number++;
           }
        }
        a2 = (Atom *)Atoms.next(i2);
      }
     a1 = (Atom *)Atoms.next(i1);
   }
   sprintf(str,"\n... %d dipolar couplings generated",number);
   log.w(str);
}