bioworld.c
Go to the documentation of this file.
1 #include "bioworld.h"
2 
3 #include "babel.h"
4 #include "link.h"
5 #include "joint.h"
6 #include "polyhedron.h"
7 #include "filename.h"
8 #include "algebra.h"
9 #include "chart.h"
10 #include "samples.h"
11 
12 #include <stdio.h>
13 #include <string.h>
14 #include <math.h>
15 
27 /****************************************************************************/
28 
49 void ReadResidueList(char *fname,unsigned int *nr,char *ch,unsigned int **r);
50 
95 void ReadRigidsAndHinges(char *fname,unsigned int **r,
96  unsigned int *nh,unsigned int **h1,unsigned int **h2,boolean **added,
97  TBioWorld *bw);
98 
111 
124 void GetAtomPositions(char *fname,double *pos,TBioWorld *bw);
125 
138 
155 void DetectLinksAndJointsFromResidues(unsigned int nr,char ch,unsigned int *rID,TBioWorld *bw);
156 
174 void DetectLinksAndJointsFromRigidsAndHinges(unsigned int *rg,
175  unsigned int nh,unsigned int *h1,unsigned int *h2,
176  TBioWorld *bw);
177 
207 void Atoms2Transforms(Tparameters *p,double *atoms,THTransform *t,
208  TBioWorld *bw);
209 
225 void InitWorldFromMolecule(Tparameters *p,double **conformation,
226  unsigned int maxAtomsLink,TBioWorld *bw);
227 
228 /********************************************************************/
229 /********************************************************************/
230 
231 void ReadResidueList(char *fname,unsigned int *nr,char *ch,unsigned int **r)
232 {
233  Tfilename fres;
234  FILE *f;
235  int out;
236  unsigned int i,m,v;
237 
238  CreateFileName(NULL,fname,NULL,RES_EXT,&fres);
239  f=fopen(GetFileFullName(&fres),"r");
240  if (!f)
241  {
242  (*nr)=0;
243  (*r)=NULL;
244  }
245  else
246  {
247  fprintf(stderr,"Reading residues from : %s\n",GetFileFullName(&fres));
248  /* first read the chain */
249  fscanf(f,"%c",ch);
250  if ((*ch>='a')&&(*ch<='z'))
251  *ch+=('A'-'a'); /* to upper */
252  if ((*ch<'A')||(*ch>'Z'))
253  Error("Wrong chain in the residue file");
254  /* now the residue identifieres */
255  (*nr)=0;
256  m=10;
257  NEW(*r,m,unsigned int);
258  do {
259  out=fscanf(f,"%u",&v);
260  if (out==1)
261  {
262  (*r)[*nr]=v;
263  (*nr)++;
264  if (*nr==m)
265  { MEM_DUP(*r,m,unsigned int); }
266  }
267  } while ((out==1)&&(out!=EOF));
268  fclose(f);
269  fprintf(stderr," Residues : %c",*ch);
270  for(i=0;i<*nr;i++)
271  fprintf(stderr," %u",(*r)[i]);
272  fprintf(stderr,"\n");
273  }
274  DeleteFileName(&fres);
275 }
276 
277 void ReadRigidsAndHinges(char *fname,unsigned int **r,
278  unsigned int *nh,unsigned int **h1,unsigned int **h2,boolean **added,
279  TBioWorld *bw)
280 {
281  Tfilename fr,fh;
282  FILE *f1,*f2;
283  int out;
284  unsigned int i,m,n,v;
285 
286  (*r)=NULL;
287  (*nh)=0;
288  (*h1)=NULL;
289  (*h2)=NULL;
290  (*added)=NULL;
291 
292  CreateFileName(NULL,fname,NULL,RIGID_EXT,&fr);
293  f1=fopen(GetFileFullName(&fr),"r");
294  CreateFileName(NULL,fname,NULL,HINGE_EXT,&fh);
295  f2=fopen(GetFileFullName(&fh),"r");
296  if ((f1)&&(f2))
297  {
298  fprintf(stderr,"Reading rigids from : %s\n",GetFileFullName(&fr));
299 
300  NEW(*r,bw->na,unsigned int);
301  /* initialy assign the atoms to no link */
302  for(i=0;i<bw->na;i++)
303  (*r)[i]=0;
304  /* scan the list of ridig assigments */
305  do {
306  out=fscanf(f1,"%u",&n); //num atom (from 1)
307  if (n>bw->na)
308  Error("Missmatch between the molecule and the rigid");
309  if (out==1)
310  {
311  out=fscanf(f1,"%u",&v); //num rigid
312  if (out==1)
313  (*r)[n-1]=v;
314  }
315  } while ((out==1)&&(out!=EOF));
316  fclose(f1);
317 
318  /* ensure that links are numbered from 1 */
319  n=0;
320  for(i=0;i<bw->na;i++)
321  {
322  if ((*r)[i]>0)
323  {
324  if ((n==0)||((*r)[i]<n))
325  n=(*r)[i];
326  }
327  }
328  if (n==0)
329  Error("Wrong rigid information in ReadRigidsAndHinges");
330  if (n>1)
331  {
332  n--;
333  /* Offset the ridig ID so that it starts at 1. */
334  for(i=0;i<bw->na;i++)
335  {
336  if ((*r)[i]>0)
337  (*r)[i]-=n;
338  }
339  }
340 
341  /* Reading the hinges */
342  fprintf(stderr,"Reading hinges from : %s\n",GetFileFullName(&fh));
343  m=10;
344  NEW(*h1,m,unsigned int);
345  NEW(*h2,m,unsigned int);
346  NEW(*added,m,boolean);
347  do {
348  out=fscanf(f2,"%u",&n); //num atom1 (from 1)
349  if (out==1)
350  {
351  out=fscanf(f2,"%u",&v); //num atoms2 (from 1)
352  if (out==1)
353  {
354  (*h1)[*nh]=n-1;
355  (*h2)[*nh]=v-1;
356 
357  /* Some of the joints might be defined over
358  unconnected atoms (h-bonds). We add them. */
359  if (!HasBond(n-1,v-1,bw->m))
360  {
361  AddBond(n-1,v-1,bw->m);
362  bw->nbAtom[n-1]++;
363  bw->nbAtom[v-1]++;
364  bw->nb+=2;
365  (*added)[*nh]=TRUE;
366  }
367  else
368  (*added)[*nh]=FALSE;
369 
370  (*nh)++;
371 
372  if (*nh==m)
373  {
374  MEM_DUP(*h1,m,unsigned int);
375  MEM_EXPAND(*h2,m,unsigned int);
376  MEM_EXPAND(*added,m,boolean);
377  }
378  }
379  }
380  } while ((out==1)&&(out!=EOF));
381  fclose(f2);
382 
383  if ((*nh)==0)
384  {
385  if ((*r)!=NULL)
386  free(*r);
387  if ((*h1)!=NULL)
388  free(*h1);
389  if ((*h2)!=NULL)
390  free(*h2);
391  (*nh)=0;
392  }
393  }
394  DeleteFileName(&fr);
395  DeleteFileName(&fh);
396 }
397 
399 {
400  unsigned int i,j;
401  TBondIterator *it;
402  boolean fix;
403 
404  bw->na=nAtoms(bw->m);
405  if (bw->na==2)
406  Error("Too small molecule (only 2 atoms!)");
407 
408  /* Number of bonds for each atom */
409  NEW(bw->nbAtom,bw->na,unsigned int);
410  for(i=0;i<bw->na;i++)
411  {
412  #if (_DEBUG>1)
413  printf("Atom %u [%g]-> ",i,VdWRadius(i,bw->m));
414  #endif
415  j=GetFirstNeighbour(i,&fix,&it,bw->m);
416  bw->nbAtom[i]=0;
417  while(j!=NO_UINT)
418  {
419  #if (_DEBUG>1)
420  if (fix)
421  printf("=%u ",j);
422  else
423  printf("-%u ",j);
424  #endif
425  bw->nbAtom[i]++;
426  j=GetNextNeighbour(i,&fix,it,bw->m);
427  }
428  DeleteBondIterator(it);
429  bw->nb+=bw->nbAtom[i];
430  bw->nba+=(bw->nbAtom[i]*(bw->nbAtom[i]-1));
431  #if (_DEBUG>1)
432  printf("\n");
433  #endif
434  }
435 }
436 
437 void GetAtomPositions(char *fname,double *pos,TBioWorld *bw)
438 {
439  Tfilename fatoms;
440  unsigned int i,j,l;
441  FILE *f;
442 
443  /* Atom positions in the pdb are given with very low accuracy
444  We try to use better estimations, if available. */
445  CreateFileName(NULL,fname,NULL,ATOM_EXT,&fatoms);
446 
447  f=fopen(GetFileFullName(&fatoms),"r");
448  if (!f)
449  GetAtomCoordinates(pos,bw->m);
450  else
451  {
452  /* the atoms file must include the 3D position of each
453  atom in the same order as listed in the pdb */
454  fprintf(stderr,"Reading atom positions from : %s\n",GetFileFullName(&fatoms));
455 
456  l=0;
457  for(i=0;i<bw->na;i++)
458  {
459  for(j=0;j<3;j++,l++)
460  fscanf(f,"%lf",&(pos[l]));
461  }
462  fclose(f);
463  }
464  DeleteFileName(&fatoms);
465 
466  #if (_DEBUG>2)
467  for(i=0,j=0;i<bw->na;i++,j+=3)
468  printf("%g %g %g\n",pos[j],pos[j+1],pos[j+2]);
469  #endif
470 }
471 
473 {
474  unsigned int i,n,k,l1,r,l;
475  boolean fix;
476  TBondIterator *it;
477 
478  NEW(bw->linkID,bw->na,unsigned int);
479  NEW(bw->linkList,bw->na,unsigned int);
480  InitVector(sizeof(unsigned int),CopyID,DeleteID,10,&(bw->links));
481  InitVector(sizeof(unsigned int),CopyID,DeleteID,10,&(bw->joint1));
482  InitVector(sizeof(unsigned int),CopyID,DeleteID,10,&(bw->joint2));
483 
484  for(i=0;i<bw->na;i++)
485  bw->linkID[i]=NO_UINT;
486 
487  bw->cut=FALSE; /* all atoms will be used to define the world */
488 
489  n=0; /* atoms assigned to link so far */
490  k=0; /* atom to process next */
491  bw->nl=0; /* Number of links so far */
492  bw->nj=0; /* Number of joints so far */
493  while(n<bw->na)
494  {
495  /* Skip atoms already assigned and isolated atoms */
496  while ((bw->linkID[k]!=NO_UINT)||(bw->nbAtom[k]<=1)) k++;
497 
498  NewVectorElement((void*)&k,&(bw->links)); /* atom 'k' starts a new link */
499  bw->linkList[k]=NO_UINT; /* mark end of list */
500  bw->linkID[k]=bw->nl;
501  n++; /* ane atom more already assigned */
502 
503  l1=k; /* last atom added to the solid */
504  r=l1; /* last expanded atom in the solid */
505 
506  /* Check if the link can be extended from any of the atoms
507  already in it -> discard */
508  while (r!=NO_UINT)
509  {
510  /* Terminal atoms are reached by other atoms in the rigid group
511  and are not a point from where to extend the rigid. */
512  if (bw->nbAtom[r]>1)
513  {
514  l=GetFirstNeighbour(r,&fix,&it,bw->m); /* 1st neigh. for atom 'r' */
515  /* Check all neighbours for atom 'r' */
516  while(l!=NO_UINT)
517  {
518  /* If not visited and the bond to the atom is double or if
519  the atom is isolated -> add to the link */
520  if ((fix)||(bw->nbAtom[l]==1))
521  {
522  if (bw->linkID[l]==NO_UINT)
523  {
524  bw->linkList[l1]=l;
525  bw->linkList[l]=NO_UINT; /* end of list */
526  bw->linkID[l]=bw->nl;
527  n++;
528 
529  l1=l; /* we have a new last atom in the link */
530  }
531  }
532  else
533  {
534  if (r<l) /* avoid defining the same joint twice */
535  {
536  NewVectorElement((void*)&r,&(bw->joint1));
537  NewVectorElement((void*)&l,&(bw->joint2));
538  bw->nj++;
539  }
540  }
541  /* move to next */
542  l=GetNextNeighbour(r,&fix,it,bw->m); /* Next neigh. for 'r' */
543  }
544  DeleteBondIterator(it);
545  }
546  /* Go for next atom in the rigid (if any) */
547  r=bw->linkList[r];
548  }
549  bw->nl++; /* we have a new link */
550  }
551 }
552 
553 void DetectLinksAndJointsFromResidues(unsigned int nr,char ch,unsigned int *rID,TBioWorld *bw)
554 {
555  unsigned int i,j,k,l1,lz,r,l,ar;
556  boolean fix,found,revolute,proline;
557  unsigned int *nID,*caID,*cID;
558  TBondIterator *it;
559  double *pos;
560 
561  if ((nr==NO_UINT)||(rID==NULL))
562  Error("Using DetectLinksAndJointsFromResidues with an empty list of residues");
563 
564  /* Ensure that the residues are increasing */
565  for(i=1;i<nr;i++)
566  {
567  if (rID[i-1]>=rID[i])
568  Error("The residue identifiers must be sorted (increasing)");
569  }
570 
571  NEW(bw->linkID,bw->na,unsigned int);
572  NEW(bw->linkList,bw->na,unsigned int);
573  InitVector(sizeof(unsigned int),CopyID,DeleteID,10,&(bw->links));
574  InitVector(sizeof(unsigned int),CopyID,DeleteID,10,&(bw->joint1));
575  InitVector(sizeof(unsigned int),CopyID,DeleteID,10,&(bw->joint2));
576 
577  for(i=0;i<bw->na;i++)
578  bw->linkID[i]=NO_UINT;
579 
580  /* Only have joints in between N-Calpha and Calpha-C fore ach residue */
581  /* Assume that the N C-alpha and C are the first 3 atoms for each residue
582  and look for them, storing their indices */
583  NEW(nID,nr,unsigned int);
584  NEW(caID,nr,unsigned int);
585  NEW(cID,nr,unsigned int);
586 
587  for(i=0;i<nr;i++)
588  {
589  nID[i]=NO_UINT;
590  caID[i]=NO_UINT;
591  cID[i]=NO_UINT;
592  }
593 
594  for(i=0;i<bw->na;i++)
595  {
596  ar=GetAtomResidue(i,bw->m);
597 
598  /* if the atom is in a residue of the correct chain */
599  if ((ar!=NO_UINT)&&(GetAtomChain(i,bw->m)==ch))
600  {
601  /* If this is one of the flexible residues ... */
602  j=0;
603  found=FALSE;
604  while((!found)&&(j<nr))
605  {
606  found=(ar==rID[j]);
607  if (!found) j++;
608  }
609 
610  if (found)
611  {
612  if (nID[j]==NO_UINT)
613  nID[j]=i; /* First atom in res. */
614  else
615  {
616  if (caID[j]==NO_UINT)
617  caID[j]=i; /* Second atom in res. */
618  else
619  {
620  if (cID[j]==NO_UINT)
621  cID[j]=i; /* Third atom in res. */
622  }
623  }
624  }
625  }
626  }
627 
628  /* Some checks to ensure that everything is fine */
629  for(i=0;i<nr;i++)
630  {
631  if ((nID[i]==NO_UINT)||(caID[i]==NO_UINT)||(cID[i]==NO_UINT))
632  Error("Undefined residue");
633  }
634 
635  /* Adjust the bounding box including the loop */
636  bw->cut=TRUE; /* only the atoms in a given box will be used to
637  define the world */
638  pos=&(bw->pos[3*nID[0]]);
639  InitBoxFromPoint(3,pos,&(bw->cutB));
640  for(i=nID[0];i<=cID[nr-1];i++)
641  {
642  pos=&(bw->pos[3*i]);
643  ExpandBox(pos,&(bw->cutB));
644  }
645  bw->cut=FALSE; /* Not actually using the bounding box .... FIX THIS */
646 
647  /* Now define the links and joints */
648 
649  /* Reserve link 0 for the non-mobile residues (link 0 is the ground link, the
650  one the defines the global reference frame).
651  Right now we just initialize it. We will add more atoms after
652  defining the mobile links. */
653  lz=nID[0]; /* This will be fist atom in link 0 */
654  NewVectorElement((void*)&lz,&(bw->links)); /* atom 'lz' starts a new link */
655  bw->linkList[lz]=NO_UINT; /* mark end of list */
656  bw->linkID[lz]=0; /* 'lz' is part of link 0 */
657  bw->nl=1; /* we have one link */
658 
659  bw->nj=0; /* Number of joints so far */
660 
661  /* Now define the mobile links, two per flexible residue, except for the last one
662  (it connects to link 0) */
663  for(i=0;i<2*nr-1;i++)
664  {
665  k=i/2; /* number in the array of flexible residues */
666  /* Select the atom from which to start the current rigid group of atoms */
667  if ((i%2)==0)
668  r=caID[k];
669  else
670  r=cID[k];
671 
672  /* Initialize the group with the selected atom */
673  NewVectorElement((void*)&r,&(bw->links)); /* atom 'r' starts a new link */
674  bw->linkList[r]=NO_UINT; /* mark end of list */
675  bw->linkID[r]=bw->nl;
676 
677  l1=r; /* last atom added to the solid */
678 
679  /* Add atoms to the group propagating over the atom bonds recursively.
680  Do not propatage through the rotable bonds. N-Ca-C */
681  while (r!=NO_UINT)
682  {
683  /* Terminal atoms are reached by other atoms in the rigid group
684  and are not a point from where to extend the rigid.
685  Isolated atoms at the end of the loop are also considered, though */
686  if (bw->nbAtom[r]>1)
687  {
688  l=GetFirstNeighbour(r,&fix,&it,bw->m); /* 1st neigh. for atom 'r' */
689 
690  /* Check all neighbours for atom 'r' */
691  while(l!=NO_UINT)
692  {
693  /* Do not propagate via the proline loop: Skip the bond closing over the
694  N atom in the backbone */
695  if ((IsAtomInProline(r,bw->m))&&(GetAtomResidue(r,bw->m)==rID[k])&&(GetAtomResidue(l,bw->m)==rID[k]))
696  proline=(((r==nID[k])&&(GetAtomicNumber(l,bw->m)==6)&&(l!=caID[k]))||
697  ((l==nID[k])&&(GetAtomicNumber(r,bw->m)==6)&&(r!=caID[k])));
698  else
699  proline=FALSE;
700 
701  if (!proline)
702  {
703  /* Propagating from a Calpha atom -> the expansion is limited by C-Ca and
704  the Ca-C bonds in the current resiude
705  Propagating from a C atom -> the expansion is limited by C-Ca bond
706  in the current residue and the N-Ca bond in the next one */
707  if (i%2==0)
708  revolute=(((r==caID[k])&&(l==nID[k]))||
709  ((r==caID[k])&&(l==cID[k])));
710  else
711  revolute=(((r==cID[k])&&(l==caID[k]))||
712  ((r==nID[k+1])&&(l==caID[k+1])));
713 
714  if (revolute)
715  {
716  /* Actually, revolute bonds define a joint */
717  if (r==caID[k]) /* but not twice! */
718  {
719  NewVectorElement((void*)&r,&(bw->joint1));
720  NewVectorElement((void*)&l,&(bw->joint2));
721  bw->nj++;
722  }
723  }
724  else
725  {
726  if (bw->linkID[l]==NO_UINT) /* atom still not assigned to any link */
727  {
728  bw->linkList[l1]=l;
729  bw->linkList[l]=NO_UINT; /* end of list */
730  bw->linkID[l]=bw->nl;
731 
732  l1=l; /* we have a new last atom in the link */
733  }
734  }
735  }
736  /* move to next bonded atom */
737  l=GetNextNeighbour(r,&fix,it,bw->m); /* Next neigh. for 'r' */
738  }
739  /* Delete the iterator over bonded atoms */
740  DeleteBondIterator(it);
741  }
742  /* Go for next atom in the rigid (if any) and propagate form them */
743  r=bw->linkList[r];
744  }
745  bw->nl++; /* we have a new link */
746  }
747 
748  /* All the non-assigned atoms must be in link 0 (even if not linked between them) */
749  for(i=0;i<bw->na;i++)
750  {
751  if (bw->linkID[i]==NO_UINT)
752  {
753  bw->linkList[lz]=i;
754  bw->linkList[i]=NO_UINT;
755  bw->linkID[i]=0;
756  lz=i;
757  }
758  }
759 
760  /* Release index of key atoms */
761  free(nID);
762  free(caID);
763  free(cID);
764 }
765 
767  unsigned int nh,unsigned int *h1,unsigned int *h2,
768  TBioWorld *bw)
769 {
770  unsigned int i,j,k,l1,lz,r,l,m,nl;
771  boolean fix,found,stop;
772  TBondIterator *it;
773 
774  if ((rg==NULL)||(nh==NO_UINT)||(nh==0)||(h1==NULL)||(h2==NULL))
775  Error("Using DetectLinksAndJointsFromRigidsAndHinges with an empty list of residues");
776 
777  NEW(bw->linkID,bw->na,unsigned int);
778  NEW(bw->linkList,bw->na,unsigned int);
779  InitVector(sizeof(unsigned int),CopyID,DeleteID,10,&(bw->links));
780  InitVector(sizeof(unsigned int),CopyID,DeleteID,10,&(bw->joint1));
781  InitVector(sizeof(unsigned int),CopyID,DeleteID,10,&(bw->joint2));
782 
783  for(i=0;i<bw->na;i++)
784  bw->linkID[i]=NO_UINT;
785 
786  bw->cut=FALSE;
787 
788  /* Now define the links and joints */
789  /* Define the given ridigs */
790  /* determine the maximum id for the given rigids */
791  m=0;
792  for(i=0;i<bw->na;i++)
793  {
794  if ((rg[i]!=NO_UINT)&&(rg[i]>m)) m=rg[i];
795  }
796  /* for all the rigids */
797  bw->nl=0;
798  for(i=1;i<=m;i++) /* rigids are num. from 1 */
799  {
800  found=FALSE;
801  lz=0;
802  while ((!found)&&(lz<bw->na))
803  {
804  found=(rg[lz]==i);
805  if (!found)
806  lz++;
807  }
808  if (found)
809  {
810  /* 'lz' is the first atom of the rigid bw->nl */
811  if (bw->linkID[lz]!=NO_UINT)
812  Error("Atom assigned to two different rigids");
813 
814  NewVectorElement((void*)&lz,&(bw->links)); /* atom 'lz' starts a new link */
815  bw->linkList[lz]=NO_UINT; /* mark end of list */
816  bw->linkID[lz]=bw->nl; /* 'lz' is part of link 0 */
817  for(j=lz+1;j<bw->na;j++)
818  {
819  if (rg[j]==i)
820  {
821  if (bw->linkID[j]!=NO_UINT)
822  Error("Atom assigned to two different rigids");
823  bw->linkList[lz]=j;
824  bw->linkList[j]=NO_UINT;
825  bw->linkID[j]=bw->nl;
826  lz=j;
827  }
828  }
829  bw->nl++;
830  }
831  }
832 
833  /* Now define the rigids from the joints. This must somehow connect
834  the given rigids */
835  nl=bw->nl;
836 
837  for(i=0;i<2*nh;i++)
838  {
839  k=i/2; /* number in the array of joints */
840  /* Select one of the atoms at the end of a joint */
841  if ((i%2)==0)
842  r=h1[k];
843  else
844  r=h2[k];
845 
846  /* if the atom is not yet in a rigid.... */
847  if (bw->linkID[r]==NO_UINT)
848  {
849  /* Initialize the group with the selected atom */
850  NewVectorElement((void*)&r,&(bw->links)); /* atom 'r' starts a new link */
851  bw->linkList[r]=NO_UINT; /* mark end of list */
852  bw->linkID[r]=bw->nl;
853 
854  l1=r; /* last atom added to the solid */
855 
856  /* Propagate from atom 'r' */
857  while (r!=NO_UINT)
858  {
859  /* Terminal atoms are reached by other atoms in the rigid group
860  and are not a point from where to extend the rigid.
861  Isolated atoms at the end of the loop are also considered, though */
862  if (bw->nbAtom[r]>1)
863  {
864  l=GetFirstNeighbour(r,&fix,&it,bw->m); /* 1st neigh. for atom 'r' */
865 
866  /* Check all neighbours for atom 'r' */
867  while(l!=NO_UINT)
868  {
869  /* stop the expansion if we reach a solid or if (l,r) is one of the
870  given joints */
871  if (bw->linkID[l]!=NO_UINT)
872  stop=TRUE;
873  else
874  {
875  stop=FALSE;
876  k=0;
877  while((!stop)&&(k<nh))
878  {
879  stop=(((h1[k]==r)&&(h2[k]==l))||((h1[k]==l)&&(h2[k]==r)));
880  k++;
881  }
882  }
883 
884  if (!stop)
885  {
886  /* atom 'l' has to be added to the link */
887  bw->linkList[l1]=l;
888  bw->linkList[l]=NO_UINT; /* end of list */
889  bw->linkID[l]=bw->nl;
890 
891  l1=l; /* we have a new last atom in the link */
892  }
893  /* move to next bonded atom */
894  l=GetNextNeighbour(r,&fix,it,bw->m); /* Next neigh. for 'r' */
895  }
896  /* Delete the iterator over bonded atoms */
897  DeleteBondIterator(it);
898  }
899  /* Go for next atom in the rigid (if any) and propagate form it */
900  r=bw->linkList[r];
901  }
902  bw->nl++; /* we have a new link */
903  }
904  }
905 
906  /* Try to expand the given rigids as much as possible to capture un-assigned atoms */
907  for(i=0;i<nl;i++)
908  {
909  r=*(unsigned int*)GetVectorElement(i,&(bw->links)); /* First atom in the link */
910  if (r==NO_UINT)
911  Error("A link without atoms?");
912 
913  /* determine the last atom added to this link */
914  j=r;
915  l1=NO_UINT;
916  while(j!=NO_UINT)
917  {
918  l1=j;
919  j=bw->linkList[j];
920  }
921 
922  /* Propagate from atom 'r' */
923  while (r!=NO_UINT)
924  {
925  if (bw->nbAtom[r]>1)
926  {
927  l=GetFirstNeighbour(r,&fix,&it,bw->m); /* 1st neigh. for atom 'r' */
928 
929  /* Check all neighbours for atom 'r' */
930  while(l!=NO_UINT)
931  {
932  /* if we reach a non-assigned atom 'l' */
933  if (bw->linkID[l]==NO_UINT)
934  {
935  /* atom 'l' has to be added to the link */
936  bw->linkList[l1]=l;
937  bw->linkList[l]=NO_UINT; /* end of list */
938  bw->linkID[l]=i;
939 
940  l1=l; /* we have a new last atom in the link */
941  }
942  /* move to next bonded atom */
943  l=GetNextNeighbour(r,&fix,it,bw->m); /* Next neigh. for 'r' */
944  }
945  /* Delete the iterator over bonded atoms */
946  DeleteBondIterator(it);
947  }
948  /* Continue propagation from the just added atoms */
949  r=bw->linkList[r];
950  }
951  }
952 
953  /* joints are given as parameters */
954  bw->nj=nh;
955  for(i=0;i<bw->nj;i++)
956  {
957  NewVectorElement((void*)&(h1[i]),&(bw->joint1));
958  NewVectorElement((void*)&(h2[i]),&(bw->joint2));
959  }
960 }
961 
962 void Atoms2Transforms(Tparameters *p,double *atoms,THTransform *t,TBioWorld *bw)
963 {
964  unsigned int l,i,j,k;
965  boolean fix;
966  TBondIterator *it;
967  double *p1,*p2,*p3;
968  double x[3],y[3],z[3];
969  boolean found;
970  /*
971  double c;
972  double v1[3],v2[3];
973  double epsilon;
974  */
975 
976  for(l=0;l<bw->nl;l++)
977  {
978  /* first atom in the link */
979  i=*(unsigned int *)GetVectorElement(l,&(bw->links));
980  j=NO_UINT;
981  k=NO_UINT;
982 
983  /* If we are in a protein try to get the pre-defined
984  reference frame */
985  if (GetAtomResidue(i,bw->m)!=NO_UINT)
986  {
987  /* Check if we are in the first rigid of the residue N-Ca-C */
988  /* We assume that atoms in a residue are given in order N-Ca-C */
989  found=FALSE;
990  while ((!found)&&(i!=NO_UINT))
991  {
992  found=((i>0)&&
993  (GetAtomicNumber(i,bw->m)==6)&& /*Ca*/
994  (GetAtomicNumber(i-1,bw->m)==7)&& /*N */
995  (GetAtomicNumber(i+1,bw->m)==6)&& /*C */
996  (GetAtomResidue(i,bw->m)==GetAtomResidue(i-1,bw->m))&&
997  (GetAtomResidue(i,bw->m)==GetAtomResidue(i+1,bw->m))&&
998  (HasBond(i,i-1,bw->m))&&
999  (HasBond(i,i+1,bw->m)));
1000  if (!found)
1001  i=bw->linkList[i];
1002  }
1003  if (found)
1004  {
1005  /* i = Ca */
1006  j=i+1; /* C */
1007  k=i-1; /* N */
1008  }
1009  else
1010  {
1011  /* Check if we are in the second ridig of the residue Ca-C-N */
1012  /* here Ca-C are in the same res. and num. consecutively but
1013  N is in next residue and non-consecutive */
1014  i=*(unsigned int *)GetVectorElement(l,&(bw->links));
1015  while ((!found)&&(i!=NO_UINT))
1016  {
1017  j=GetFirstNeighbour(i,&fix,&it,bw->m);
1018  while ((!found)&&(j!=NO_UINT))
1019  {
1020  found=((i>0)&&
1021  (GetAtomicNumber(i,bw->m)==6)&& /*C */
1022  (GetAtomicNumber(i-1,bw->m)==6)&& /*Ca*/
1023  (GetAtomicNumber(j,bw->m)==7)&& /*N */
1024  (GetAtomResidue(i,bw->m)==GetAtomResidue(i-1,bw->m))&&
1025  (GetAtomResidue(i,bw->m)+1==GetAtomResidue(j,bw->m))&&
1026  (HasBond(i,i-1,bw->m))&&
1027  (HasBond(i,j,bw->m)));
1028  if (!found)
1029  j=GetNextNeighbour(i,&fix,it,bw->m);
1030  }
1031  DeleteBondIterator(it);
1032  if (!found)
1033  i=bw->linkList[i];
1034  }
1035  if (found)
1036  {
1037  /* i = Ca */
1038  /* j = N */
1039  k=i-1; /* C*/
1040  }
1041  }
1042  }
1043 
1044  /* if the previous strategy did not work just pick any 3 atoms
1045  in the link */
1046  if ((i==NO_UINT)||(j==NO_UINT)||(k==NO_UINT))
1047  {
1048  i=*(unsigned int *)GetVectorElement(l,&(bw->links));
1049  while((i!=NO_UINT)&&(bw->nbAtom[i]<2))
1050  i=bw->linkList[i];
1051 
1052  if (i==NO_UINT)
1053  Error("A link with no valid atom to define a reference frame");
1054 
1055  /* Note that 'j' and 'k' are not necessarily in the same link (=rigid
1056  group of atoms) as atom 'i'. However, the bond between them define
1057  a rigid angle and, thus, a fixed plane with respect to the link. */
1058 
1059  /* We first try to avoid using isolated atoms as references */
1060  j=GetFirstNeighbour(i,&fix,&it,bw->m);
1061  if (j==NO_UINT)
1062  Error("Inconsistency in BioWorld (Atoms2Transforms)");
1063  while ((j!=NO_UINT)&&(bw->nbAtom[j]<2))
1064  j=GetNextNeighbour(i,&fix,it,bw->m);
1065  if (j!=NO_UINT)
1066  {
1067  k=GetNextNeighbour(i,&fix,it,bw->m);
1068  while ((k!=NO_UINT)&&(bw->nbAtom[k]<2))
1069  k=GetNextNeighbour(i,&fix,it,bw->m);
1070  }
1071  DeleteBondIterator(it);
1072 
1073  /* If not possible just use the first two atoms */
1074  if ((j==NO_UINT)||(k==NO_UINT))
1075  {
1076  j=GetFirstNeighbour(i,&fix,&it,bw->m);
1077  k=GetNextNeighbour(i,&fix,it,bw->m);
1078  if (k==NO_UINT)
1079  Error("Inconsistency in BioWorld (Atoms2Transforms)");
1080  DeleteBondIterator(it);
1081  }
1082  }
1083 
1084  p1=&(atoms[3*i]);
1085  p2=&(atoms[3*j]);
1086  p3=&(atoms[3*k]);
1087 
1088  DifferenceVector(3,p2,p1,x);
1089  Normalize(3,x);
1090 
1091  DifferenceVector(3,p3,p1,y);
1092  Normalize(3,y);
1093 
1094  CrossProduct(x,y,z);
1095  Normalize(3,z);
1096 
1097  CrossProduct(z,x,y);
1098  Normalize(3,y); /* just for numerical accuracy */
1099 
1100  /*
1101  epsilon=GetParameter(CT_EPSILON,p);
1102 
1103  DifferenceVector(3,p2,p1,v1);
1104  Normalize(3,v1);
1105 
1106  DifferenceVector(3,p3,p1,v2);
1107  Normalize(3,v2);
1108 
1109  memcpy(x,v1,3*sizeof(double));
1110 
1111  c=DotProduct(v1,v2);
1112  if (fabs(1-fabs(c))<epsilon)
1113  Error("Aligned bonds");
1114  SumVectorScale(3,v2,-c,v1,y);
1115  ScaleVector(1.0/sqrt(1-c*c),3,y);
1116 
1117  CrossProduct(x,y,z);
1118  */
1119 
1120  HTransformFromVectors(x,y,z,p1,&(t[l]));
1121 
1122 #if (0)
1123  fprintf(stderr,"Link %u: [%u %u %u]\n",l,i,j,k);
1124  if (l==0)
1125  HTransformPrint(stderr,&(t[l]));
1126 #endif
1127  }
1128 }
1129 
1130 void InitWorldFromMolecule(Tparameters *p,double **conformation,
1131  unsigned int maxAtomsLink,TBioWorld *bw)
1132 {
1133  unsigned int rep,i,j,m,n;
1134  double *x,*y,*z;
1135  boolean fix,simple;
1136  Tpolyhedron sphere,cylinder,segments;
1137  Tlink link;
1138  Tjoint joint;
1139  unsigned int k,l,nal;
1140  Tcolor defaultColor,carbon,sulphur,hydrogen,black,oxygen,nitrogen,red,green,blue,*atomColor;
1141  char name[100];
1142  double **jointPoint;
1143  double vdwRatio; /* ratio of Van der Waals radius used */
1144  double radius,*atomPos,endPoint[3];
1145  TBondIterator *it;
1146  THTransform *l2g;
1147 
1148  for(i=0;i<bw->na;i++)
1149  {
1150  if (bw->linkID[i]==NO_UINT)
1151  Error("Unassigned atom to link in InitWorldFromMolecule (this atom can not be correctly moved)");
1152  }
1153 
1154  /* Get the transform to global coordinates for each link. */
1155  NEW(l2g,bw->nl,THTransform);
1156  Atoms2Transforms(p,bw->pos,l2g,bw);
1157 
1158  /* Compute the global2local (inverse of local2global) */
1159  NEW(bw->g2l,bw->nl,THTransform);
1160  for(i=0;i<bw->nl;i++)
1161  HTransformInverse(&(l2g[i]),&(bw->g2l[i]));
1162 
1163  /* Use the transforms to compute the local coordinates of each atom */
1164  NEW(bw->localPos,bw->na*3,double);
1165  for(i=0,l=0;i<bw->na;i++,l+=3)
1166  {
1167  if (bw->linkID[i]!=NO_UINT) /* Only for actually used atoms */
1168  HTransformApply(&(bw->pos[l]),&(bw->localPos[l]),&(bw->g2l[bw->linkID[i]]));
1169  }
1170 
1171  /*Generate an empty world*/
1172  InitWorld(&(bw->w));
1173 
1174  /* Add the links to the world */
1175  NewColor(0.75,0.75,0.75,&defaultColor);
1176  NewColor(0.25,0.25,0.25,&carbon);
1177  NewColor(0.6,0.6,0.2,&sulphur);
1178  NewColor(1,1,1,&hydrogen);
1179  NewColor(1,0,0,&oxygen);
1180  NewColor(0,0,1,&nitrogen);
1181  NewColor(1,0,0,&red);
1182  NewColor(0,1,0,&green);
1183  NewColor(0,0,1,&blue);
1184  NewColor(0,0,0,&black);
1185  vdwRatio=GetParameter(CT_VDW_RATIO,p);
1186 
1187  for(i=0;i<bw->nl;i++)
1188  {
1189  /* Initialize the link */
1190  sprintf(name,"link_%u",i);
1191  InitLink(name,&link);
1192 
1193  /* Define a sphere for each atom in the link */
1194  /* Count the num. atoms in the link. */
1195  k=*(unsigned int *)GetVectorElement(i,&(bw->links));
1196  nal=0; /* number atoms link */
1197  while (k!=NO_UINT)
1198  {
1199  nal++;
1200  k=bw->linkList[k];
1201  }
1202  /* If we have many atoms we use a simple representation */
1203  simple=((maxAtomsLink!=NO_UINT)&&(nal>=maxAtomsLink));
1204 
1205  if (simple)
1206  {
1207  n=0;
1208  m=100;
1209  NEW(x,m,double);
1210  NEW(y,m,double);
1211  NEW(z,m,double);
1212  }
1213 
1214  /* first atom in the link */
1215  k=*(unsigned int *)GetVectorElement(i,&(bw->links));
1216  do{
1217  /* Only consider atoms in a range +/-10 atoms around the cut
1218  box, if any. */
1219  if ((!bw->cut)||(PointInBox(NULL,3,&(bw->pos[3*k]),GetBoxDiagonal(NULL,&(bw->cutB))/2,&(bw->cutB))))
1220  {
1221  atomPos=&(bw->localPos[3*k]);
1222 
1223  radius=vdwRatio*VdWRadius(k,bw->m);
1224 
1225  switch(GetAtomicNumber(k,bw->m))
1226  {
1227  case 1:
1228  atomColor=&hydrogen;
1229  break;
1230  case 6:
1231  atomColor=&carbon;
1232  break;
1233  case 7:
1234  atomColor=&nitrogen;
1235  break;
1236  case 8:
1237  atomColor=&oxygen;
1238  break;
1239  case 16:
1240  atomColor=&sulphur;
1241  break;
1242  default:
1243  atomColor=&defaultColor;
1244  }
1245  if (simple)
1246  NewSphere(radius,atomPos,atomColor,0,HIDDEN_SHAPE,&sphere);
1247  else
1248  NewSphere(radius,atomPos,atomColor,0,NORMAL_SHAPE,&sphere);
1249  AddBody2Link(&sphere,&link);
1250  DeletePolyhedron(&sphere);
1251 
1252  /* A cylinder/line for each pair of bonded atoms. */
1253  j=GetFirstNeighbour(k,&fix,&it,bw->m);
1254  while(j!=NO_UINT)
1255  {
1256  if (simple)
1257  {
1258  if ((k<j)||(bw->linkID[j]!=bw->linkID[k]))
1259  {
1260  /* the end point might be in another link and we
1261  need its coordinates in this link */
1262  HTransformApply(&(bw->pos[3*j]),endPoint,&(bw->g2l[i]));
1263 
1264  if (n>=m-1)
1265  {
1266  MEM_DUP(x,m,double);
1267  MEM_EXPAND(y,m,double);
1268  MEM_EXPAND(z,m,double);
1269  }
1270  x[n]=atomPos[0];
1271  y[n]=atomPos[1];
1272  z[n]=atomPos[2];
1273  n++;
1274  x[n]=endPoint[0];
1275  y[n]=endPoint[1];
1276  z[n]=endPoint[2];
1277  n++;
1278  }
1279  }
1280  else
1281  {
1282  /* the end point might be in another link and we
1283  need its coordinates in this link */
1284  HTransformApply(&(bw->pos[3*j]),endPoint,&(bw->g2l[i]));
1285  DifferenceVector(3,endPoint,atomPos,endPoint);
1286  SumVectorScale(3,atomPos,0.5,endPoint,endPoint);
1287 
1288  NewCylinder(0.15,atomPos,endPoint,atomColor,0,DECOR_SHAPE,&cylinder);
1289  AddBody2Link(&cylinder,&link);
1290  DeletePolyhedron(&cylinder);
1291  }
1292  j=GetNextNeighbour(k,&fix,it,bw->m);
1293  }
1294  DeleteBondIterator(it);
1295  }
1296 
1297  /* Go for next atom */
1298  k=bw->linkList[k];
1299  } while(k!=NO_UINT);
1300 
1301  if (simple)
1302  {
1303  switch(i%4)
1304  {
1305  case 0:
1306  NewSegments(n,x,y,z,&blue,&segments);
1307  break;
1308  case 1:
1309  NewSegments(n,x,y,z,&green,&segments);
1310  break;
1311  case 2:
1312  NewSegments(n,x,y,z,&red,&segments);
1313  break;
1314  case 3:
1315  NewSegments(n,x,y,z,&black,&segments);
1316  break;
1317  }
1318  AddBody2Link(&segments,&link);
1319  DeletePolyhedron(&segments);
1320  free(x);
1321  free(y);
1322  free(z);
1323  }
1324 
1325  AddLink2World(&link,FALSE,&(bw->w));
1326  DeleteLink(&link);
1327  }
1328  DeleteColor(&defaultColor);
1329  DeleteColor(&carbon);
1330  DeleteColor(&hydrogen);
1331  DeleteColor(&sulphur);
1332  DeleteColor(&black);
1333 
1334  /* Add the joints to the world */
1335  NEW(jointPoint,4,double*);
1336  for(i=0;i<4;i++)
1337  NEW(jointPoint[i],3,double);
1338 
1339  rep=(unsigned int)(GetParameter(CT_REPRESENTATION,p));
1340 
1341  for(i=0;i<bw->nj;i++)
1342  {
1343  k=*(unsigned int *)GetVectorElement(i,&(bw->joint1));
1344  l=*(unsigned int *)GetVectorElement(i,&(bw->joint2));
1345 
1346  HTransformApply(&(bw->pos[3*k]),jointPoint[0],
1347  &(bw->g2l[bw->linkID[k]]));
1348  HTransformApply(&(bw->pos[3*l]),jointPoint[1],
1349  &(bw->g2l[bw->linkID[k]]));
1350 
1351  HTransformApply(&(bw->pos[3*k]),jointPoint[2],
1352  &(bw->g2l[bw->linkID[l]]));
1353  HTransformApply(&(bw->pos[3*l]),jointPoint[3],
1354  &(bw->g2l[bw->linkID[l]]));
1355 
1356  NewRevoluteJoint(i,rep,
1357  bw->linkID[k],GetWorldLink(bw->linkID[k],&(bw->w)),
1358  bw->linkID[l],GetWorldLink(bw->linkID[l],&(bw->w)),
1359  jointPoint,FALSE,NULL,NULL,FALSE,0,NULL,
1360  &joint);
1361 
1362  AddJoint2World(&joint,&(bw->w));
1363  DeleteJoint(&joint);
1364  }
1365 
1366  /* Release temporal information used in the world definition. */
1367  for(i=0;i<4;i++)
1368  free(jointPoint[i]);
1369  free(jointPoint);
1370 
1371  /* Post-process the world */
1372  if (bw->nl>0)
1373  {
1374  /* Check all collisions (between atoms in different rigids/links) */
1375  CheckAllCollisions(0,0,&(bw->w));
1376  /* Generate the kinematic equations */
1377  GenerateWorldEquations(p,&(bw->w));
1378  }
1379 
1380  /* Now that the mechanism is alraedy defined, the local2global
1381  transforms can be used to deduce the conformation in
1382  internal coordinates*/
1383  GetSolutionPointFromLinkTransforms(p,l2g,NULL,conformation,&(bw->w));
1384 
1385  /* Release the local 2 global transforms */
1386  for(i=0;i<bw->nl;i++)
1387  HTransformDelete(&(l2g[i]));
1388  free(l2g);
1389 }
1390 
1391 /********************************************************************/
1392 /********************************************************************/
1393 
1394 boolean InitBioWorld(Tparameters *p,char *filename,unsigned int maxAtomsLink,
1395  double **conformation,TBioWorld *bw)
1396 {
1397 
1398  unsigned int nr,*r;
1399  char ch;
1400  double *c;
1401  unsigned int nh,*h1,*h2;
1402  boolean *added;
1403 
1404  /* Read the molecule */
1405  bw->m=ReadMolecule(filename);
1406  if (bw->m!=NULL)
1407  {
1408  fprintf(stderr,"Reading bio-info from : %s\n",filename);
1409 
1410  /* Collect information about the molecule and eventually
1411  print debug information about it */
1413 
1414  /* Atom positions */
1415  NEW(bw->pos,3*bw->na,double);
1416  GetAtomPositions(filename,bw->pos,bw);
1417  SetAtomCoordinates(bw->pos,bw->m);
1418 
1419  /* Detect the links and joints
1420  links=rigid groups of atoms
1421  joints=revolute joints along single bonds
1422  */
1423  ReadRigidsAndHinges(filename,&r,&nh,&h1,&h2,&added,bw);
1424  if ((r!=NULL)&&(h1!=NULL)&&(h1!=NULL)&&(added!=NULL))
1426  else
1427  {
1428  ReadResidueList(filename,&nr,&ch,&r);
1429  if (r!=NULL)
1430  {
1432  free(r);
1433  }
1434  else
1436  }
1437 
1438  /* Define the world from the molecular information defined so far */
1439  InitWorldFromMolecule(p,conformation,maxAtomsLink,bw);
1440 
1441  //AdjustBioWorldGeometry(p,bw);
1442 
1443  /* Set the atom positions in accordance to the new reference frames
1444  just defined */
1447  free(c);
1448 
1449  if ((r!=NULL)&&(h1!=NULL)&&(h1!=NULL)&&(added!=NULL))
1450  {
1451  unsigned int i;
1452 
1453  /* Remove the artificially created bonds, if any. */
1454  for (i=0;i<nh;i++)
1455  {
1456  if (added[i])
1457  {
1458  RemoveBond(h1[i],h2[i],bw->m);
1459  bw->nbAtom[h1[i]]--;
1460  bw->nbAtom[h2[i]]--;
1461  bw->nb-=2;
1462  }
1463  }
1464  free(r);
1465  free(h1);
1466  free(h2);
1467  free(added);
1468  }
1469  return(TRUE);
1470  }
1471  else
1472  return(FALSE);
1473 }
1474 
1476 {
1477  TCuikSystem c;
1478  Tvariable var;
1479  Tequation eq,eq2;
1480  Tmonomial m,m2;
1481  double *s;
1482  char *vname;
1483  unsigned int i,j,l,n,nv,k;
1484  boolean fix;
1485  Tinterval range;
1486  unsigned int vID[3],vID2[3],lID;
1487  unsigned int an1,an2,an3;
1488  TBondIterator *it,*it2;
1489  double v1[3],v2[3];
1490  double sgn;
1491 
1492  InitCuikSystem(&c);
1493 
1494  nv=bw->na+bw->nb+bw->nba;
1495  NEW(s,nv,double);
1496  k=0;
1497 
1498  /* x,y,z for each atom */
1499  NEW(vname,100,char);
1500  InitMonomial(&m);
1501  InitMonomial(&m2);
1502 
1503  for(i=0;i<bw->na;i++)
1504  {
1505  for(j=0;j<3;j++)
1506  {
1507  /* The ct value from the bio-info (initial approx.
1508  to the solution) */
1509  s[k]=bw->pos[k];
1510 
1511  /* x,y,z position for each atom */
1512  sprintf(vname,"p_%u_%s",i,(j==0?"x":(j==1?"y":"z")));
1513 
1514  /* add the new variable to the system */
1515  NewVariable(SYSTEM_VAR,vname,&var);
1516  //NewInterval(s[k]-0.1,s[k]+0.1,&range);
1517  NewInterval(-20,20,&range);
1518  SetVariableInterval(&range,&var);
1519  lID=AddVariable2CS(&var,&c);
1520  DeleteVariable(&var);
1521 
1522  /* fix first atom */
1523  if (i==0)
1524  {
1525  InitEquation(&eq);
1527  SetEquationCmp(EQU,&eq);
1528  //SetEquationValue(s[k],&eq);
1529  SetEquationValue(0,&eq);
1530 
1531  AddVariable2Monomial(NFUN,lID,1,&m);
1532  AddMonomial(&m,&eq);
1533  ResetMonomial(&m);
1534 
1535  AddEquation2CS(p,&eq,&c);
1536  DeleteEquation(&eq);
1537  }
1538  k++;
1539  }
1540  }
1541 
1542  /* vector for each pair of bonded atoms */
1543  for(i=0;i<bw->na;i++)
1544  {
1545  l=GetFirstNeighbour(i,&fix,&it,bw->m);
1546  while(l!=NO_UINT)
1547  {
1548  if (l>i)
1549  {
1550  for(j=0;j<3;j++)
1551  {
1552  /* The ct value from the atom positions */
1553  s[k]=bw->pos[3*l+j]-bw->pos[3*i+j];
1554 
1555  /* v_i_l = p_l-p_i */
1556  sprintf(vname,"v_%u_%u_%s",i,l,(j==0?"x":(j==1?"y":"z")));
1557 
1558  /* Add the new variable v_i_l to the system */
1559  NewVariable(DUMMY_VAR,vname,&var);
1560  //NewInterval(s[k]-0.2,s[k]+0.2,&range);
1561  NewInterval(-20,20,&range);
1562  SetVariableInterval(&range,&var);
1563  vID[j]=AddVariable2CS(&var,&c);
1564  DeleteVariable(&var);
1565 
1566  /* Initialize the equation */
1567  InitEquation(&eq);
1569  SetEquationCmp(EQU,&eq);
1570  SetEquationValue(0,&eq);
1571 
1572  /* p_l */
1573  AddVariable2Monomial(NFUN,3*l+j,1,&m);
1574  AddMonomial(&m,&eq);
1575  ResetMonomial(&m);
1576 
1577  /* -p_i */
1578  AddCt2Monomial(-1.0,&m);
1579  AddVariable2Monomial(NFUN,3*i+j,1,&m);
1580  AddMonomial(&m,&eq);
1581  ResetMonomial(&m);
1582 
1583  /* -v_i_l */
1584  AddCt2Monomial(-1.0,&m);
1585  AddVariable2Monomial(NFUN,vID[j],1,&m);
1586  AddMonomial(&m,&eq);
1587  ResetMonomial(&m);
1588 
1589  AddEquation2CS(p,&eq,&c);
1590  DeleteEquation(&eq);
1591 
1592  k++;
1593  }
1594 
1595  /* The norm of v_i_l must the the same for all bonds
1596  of this type (type given by the atom types) */
1597  GenerateNormEquation(vID[0],vID[1],vID[2],0.0,&eq);
1598 
1599  an1=GetAtomicNumber(i,bw->m);
1600  an2=GetAtomicNumber(l,bw->m);
1601  if (an1<an2)
1602  sprintf(vname,"bl_%u_%u",an1,an2);
1603  else
1604  sprintf(vname,"bl_%u_%u",an2,an1);
1605  lID=GetCSVariableID(vname,&c);
1606  if (lID==NO_UINT)
1607  {
1608  /* and define the initial estimation from the atom pos */
1609  s[k]=Distance(3,&(bw->pos[3*i]),&(bw->pos[3*l]));
1610  s[k]*=s[k];
1611 
1612  /* First bond of this type found -> add the new variable */
1613  NewVariable(DUMMY_VAR,vname,&var);
1614  //NewInterval(s[k]-0.1,s[k]+0.1,&range);
1615  NewInterval(-400,400,&range);
1616  SetVariableInterval(&range,&var);
1617  lID=AddVariable2CS(&var,&c);
1618  DeleteVariable(&var);
1619 
1620  InitEquation(&eq2);
1621  SetEquationType(SYSTEM_EQ,&eq2);
1622  SetEquationCmp(EQU,&eq2);
1623  SetEquationValue(s[k],&eq2);
1624 
1625  AddVariable2Monomial(NFUN,lID,1,&m2);
1626  AddMonomial(&m2,&eq2);
1627  ResetMonomial(&m2);
1628 
1629  AddEquation2CS(p,&eq2,&c);
1630  DeleteEquation(&eq2);
1631 
1632  k++;
1633  }
1634 
1635  AddCt2Monomial(-1.0,&m);
1636  AddVariable2Monomial(NFUN,lID,1,&m);
1637  AddMonomial(&m,&eq);
1638  ResetMonomial(&m);
1639 
1640  AddEquation2CS(p,&eq,&c);
1641  DeleteEquation(&eq);
1642  }
1643  l=GetNextNeighbour(i,&fix,it,bw->m);
1644  }
1645  DeleteBondIterator(it);
1646  }
1647 
1648  /* angle between for each triplet of bonded atoms.
1649  In particular we look for all angles with vertex at atom 'i'
1650  (avoiding repetitions). */
1651  for(i=0;i<bw->na;i++)
1652  {
1653  l=GetFirstNeighbour(i,&fix,&it,bw->m);
1654  while(l!=NO_UINT)
1655  {
1656  n=GetFirstNeighbour(i,&fix,&it2,bw->m);
1657  while(n!=NO_UINT)
1658  {
1659  if (l>n)
1660  {
1661  /* Equation v1*v2'-n1*n2*cos(alpha)=0 */
1662  an1=GetAtomicNumber(i,bw->m);
1663  an2=GetAtomicNumber(l,bw->m);
1664  an3=GetAtomicNumber(n,bw->m);
1665 
1666  /* v1*v2 */
1667  sgn=1.0;
1668  for(j=0;j<3;j++)
1669  {
1670  if (l>i)
1671  sprintf(vname,"v_%u_%u_%s",i,l,(j==0?"x":(j==1?"y":"z")));
1672  else
1673  {
1674  sprintf(vname,"v_%u_%u_%s",l,i,(j==0?"x":(j==1?"y":"z")));
1675  if (j==0) sgn*=-1.0;
1676  }
1677  vID[j]=GetCSVariableID(vname,&c);
1678  if (vID[j]==NO_UINT)
1679  Error("Unknow atom difference when fixing an angle bond");
1680 
1681  if (n>i)
1682  sprintf(vname,"v_%u_%u_%s",i,n,(j==0?"x":(j==1?"y":"z")));
1683  else
1684  {
1685  sprintf(vname,"v_%u_%u_%s",n,i,(j==0?"x":(j==1?"y":"z")));
1686  if (j==0) sgn*=-1.0;
1687  }
1688  vID2[j]=GetCSVariableID(vname,&c);
1689  if (vID2[j]==NO_UINT)
1690  Error("Unknow atom difference when fixing an angle bond");
1691  }
1692 
1693  /* This initializes the equation */
1694  GenerateDotProductEquation(vID[0],vID[1],vID[2],
1695  vID2[0],vID2[1],vID2[2],
1696  NO_UINT,0,&eq);
1697 
1698  if ((an1<an2)&&(an1<an3))
1699  {
1700  /* an1 smaller */
1701  if (an2<an3)
1702  sprintf(vname,"ba_%u_%u_%u",an1,an2,an3);
1703  else
1704  sprintf(vname,"ba_%u_%u_%u",an1,an3,an2);
1705  }
1706  else
1707  {
1708  /* an2 smaller */
1709  if ((an2<an1)&&(an2<an3))
1710  {
1711  if (an1<an3)
1712  sprintf(vname,"ba_%u_%u_%u",an2,an1,an3);
1713  else
1714  sprintf(vname,"ba_%u_%u_%u",an2,an3,an1);
1715  }
1716  else
1717  {
1718  /* an3 smaller */
1719  if (an1<an2)
1720  sprintf(vname,"ba_%u_%u_%u",an3,an1,an2);
1721  else
1722  sprintf(vname,"ba_%u_%u_%u",an3,an2,an1);
1723  }
1724  }
1725 
1726  lID=GetCSVariableID(vname,&c);
1727  if (lID==NO_UINT)
1728  {
1729  /* and define the initial estimation from the current atom positions */
1730  DifferenceVector(3,&(bw->pos[3*l]),&(bw->pos[3*i]),v1);
1731  DifferenceVector(3,&(bw->pos[3*n]),&(bw->pos[3*i]),v2);
1732  s[k]=DotProduct(v1,v2);
1733 
1734  /* First angle of this type found */
1735  NewVariable(SYSTEM_VAR,vname,&var);
1736  //NewInterval(s[k]-0.1,s[k]+0.1,&range);
1737  NewInterval(-1,1,&range);
1738  SetVariableInterval(&range,&var);
1739  lID=AddVariable2CS(&var,&c);
1740  DeleteVariable(&var);
1741 
1742  InitEquation(&eq2);
1743  SetEquationType(DUMMY_EQ,&eq2);
1744  SetEquationCmp(EQU,&eq2);
1745  SetEquationValue(s[k],&eq2);
1746 
1747  AddVariable2Monomial(NFUN,lID,1,&m2);
1748  AddMonomial(&m2,&eq2);
1749  ResetMonomial(&m2);
1750 
1751  AddEquation2CS(p,&eq2,&c);
1752  DeleteEquation(&eq2);
1753 
1754  k++;
1755  }
1756 
1757  AddCt2Monomial(-sgn,&m);
1758  AddVariable2Monomial(NFUN,lID,1,&m);
1759  AddMonomial(&m,&eq);
1760  ResetMonomial(&m);
1761 
1762  AddEquation2CS(p,&eq,&c);
1763  DeleteEquation(&eq);
1764  }
1765  n=GetNextNeighbour(i,&fix,it2,bw->m);
1766  }
1767  DeleteBondIterator(it2);
1768  l=GetNextNeighbour(i,&fix,it,bw->m);
1769  }
1770  DeleteBondIterator(it);
1771  }
1772  DeleteMonomial(&m);
1773 
1774  /* Print Cuiksystem?? */
1775  {
1776  FILE *fout;
1777  Tbox b;
1778  char **varNames;
1779  boolean *systemVars;
1780 
1781  fout=fopen("geometry.cuik","w");
1782  if (!fout)
1783  Error("Could not open the file to store the geometry equations.");
1784  PrintCuikSystem(p,fout,&c);
1785  fclose(fout);
1786 
1787  InitBoxFromPoint(k,s,&b);
1788  fout=fopen("geometry.sol","w");
1789  if (!fout)
1790  Error("Could not open the file to store the initial approx.");
1791  NEW(varNames,k,char*);
1792  GetCSVariableNames(varNames,&c);
1793  GetCSSystemVars(&systemVars,&c);
1794  PrintBoxSubset(fout,systemVars,varNames,&b);
1795  free(varNames);
1796  free(systemVars);
1797  fclose(fout);
1798  DeleteBox(&b);
1799  }
1800  /* CuikNewton on Cuiksystem? */
1801 
1802  free(s);
1803  free(vname);
1804  DeleteCuikSystem(&c);
1805 }
1806 
1808 {
1809  return(&(bw->w));
1810 }
1811 
1812 unsigned int BioWorldNAtoms(TBioWorld *bw)
1813 {
1814  return(bw->na);
1815 }
1816 
1818 {
1819  return(GetWorldNumSystemVariables(&(bw->w)));
1820 }
1821 
1822 void BioWordGetAtomPositionsFromConformation(Tparameters *p,boolean simp,double *conformation,
1823  double *pos,TBioWorld *bw)
1824 {
1825  THTransform *tl,*t;
1826  TLinkConf *def;
1827  unsigned int i,k;
1828 
1829  /* Compute the position of the atoms form the conformation */
1830  GetLinkTransformsFromSolutionPoint(p,simp,conformation,&tl,&def,&(bw->w));
1831 
1832  /* for all links */
1833  for(i=0;i<bw->nl;i++)
1834  {
1835  /* Transform for this link */
1836  t=&(tl[i]);
1837  /* apply it for all atoms in this link */
1838  k=*(unsigned int *)GetVectorElement(i,&(bw->links));
1839  do{
1840  /* move the atom */
1841  HTransformApply(&(bw->localPos[3*k]),&(pos[3*k]),t);
1842 
1843  /* go for next atom in this link */
1844  k=bw->linkList[k];
1845 
1846  } while(k!=NO_UINT);
1847  }
1848 
1849  DeleteLinkTransforms(tl,def,&(bw->w));
1850 }
1851 
1852 void BioWordSetAtomPositionsFromConformation(Tparameters *p,boolean simp,double *conformation,
1853  TBioWorld *bw)
1854 {
1855  BioWordGetAtomPositionsFromConformation(p,simp,conformation,bw->pos,bw);
1856 
1857  SetAtomCoordinates(bw->pos,bw->m);
1858 }
1859 
1860 double BioWorldRMSE(double *pos,TBioWorld *bw)
1861 {
1862  double e;
1863  unsigned int i,n,k;
1864 
1865  e=0.0;
1866  n=0;
1867  for(i=0,k=0;i<bw->na;i++,k+=3)
1868  {
1869  /* atoms in link 0 do not move. When modelling protein loops
1870  link 0 can be quite large. */
1871  if (bw->linkID[i]!=0)
1872  {
1873  e+=Distance(3,&(pos[k]),&(bw->pos[k]));
1874  n++;
1875  }
1876  }
1877 
1878  return(sqrt(e/n));
1879 }
1880 
1881 
1883  double **conformation,TBioWorld *bw)
1884 {
1885  THTransform *l2g;
1886  unsigned int i,nv;
1887 
1888  /* Get the transform to global coordinates for each link. */
1889  NEW(l2g,bw->nl,THTransform);
1890  Atoms2Transforms(p,atoms,l2g,bw);
1891 
1892  nv=GetSolutionPointFromLinkTransforms(p,l2g,NULL,conformation,&(bw->w));
1893 
1894  /* Release the local-2-global transforms */
1895  for(i=0;i<bw->nl;i++)
1896  HTransformDelete(&(l2g[i]));
1897  free(l2g);
1898 
1899  return(nv);
1900 }
1901 
1902 double BioWorldEnergy(Tparameters *p,boolean simp,double *conformation,void *bw)
1903 {
1904 
1905  BioWordSetAtomPositionsFromConformation(p,simp,conformation,(TBioWorld *)bw);
1906 
1907  /* And compute the energy */
1908  return(ComputeEnergy(((TBioWorld *)bw)->m));
1909 }
1910 
1911 void SaveBioWorldBioInfo(Tparameters *p,char *fname,boolean simp,double *conformation,TBioWorld *bw)
1912 {
1913  BioWordSetAtomPositionsFromConformation(p,simp,conformation,bw);
1914 
1915  /* and save in bio-format*/
1916  WriteMolecule(fname,bw->m);
1917 }
1918 
1919 void PrintBioWorld(Tparameters *p,char *fname,int argc,char **arg,TBioWorld *bw)
1920 {
1921  PrintWorld(fname,argc,arg,&(bw->w));
1922 }
1923 
1925 {
1926  unsigned int i;
1927 
1928  DeleteWorld(&(bw->w));
1929  DeleteMolecule(bw->m);
1930 
1931  free(bw->localPos);
1932  free(bw->pos);
1933 
1934  free(bw->nbAtom);
1935 
1936  for(i=0;i<bw->nl;i++)
1937  HTransformDelete(&(bw->g2l[i]));
1938  free(bw->g2l);
1939 
1940  free(bw->linkList);
1941  DeleteVector((void *)&(bw->links));
1942  DeleteVector((void *)&(bw->joint1));
1943  DeleteVector((void *)&(bw->joint2));
1944  free(bw->linkID);
1945 
1946  if (bw->cut)
1947  DeleteBox(&(bw->cutB));
1948 }
1949 
unsigned int * linkList
Definition: bioworld.h:45
boolean HasBond(unsigned int na1, unsigned int na2, TMolecule *ml)
Determines is a given bond exists.
Definition: babel.cpp:184
boolean PointInBox(boolean *used, unsigned int n, double *v, double tol, Tbox *b)
Checks if a point is included in a(sub-) box.
Definition: box.c:356
void DeleteMolecule(TMolecule *ml)
Destructor.
Definition: babel.cpp:353
unsigned int nba
Definition: bioworld.h:33
void DeleteVector(void *vector)
Destructor.
Definition: vector.c:389
void SetAtomCoordinates(double *pos, TMolecule *ml)
Sets the positions of the atoms in the molecule.
Definition: babel.cpp:325
#define SYSTEM_EQ
One of the possible type of equations.
Definition: equation.h:146
unsigned int nAtoms(TMolecule *ml)
Number of atoms in a molecule.
Definition: babel.cpp:100
void NewSegments(unsigned int n, double *x, double *y, double *z, Tcolor *c, Tpolyhedron *p)
Constructor.
Definition: polyhedron.c:1098
void DetectLinksAndJoints(TBioWorld *bw)
Determines the rigid groups of atoms and the connections between them.
Definition: bioworld.c:472
#define FALSE
FALSE.
Definition: boolean.h:30
void DeleteBioWorld(TBioWorld *bw)
Destructor.
Definition: bioworld.c:1924
double BioWorldEnergy(Tparameters *p, boolean simp, double *conformation, void *bw)
Computes the energy of a given configuration.
Definition: bioworld.c:1902
void HTransformApply(double *p_in, double *p_out, THTransform *t)
Multiply a homogeneous transform and a vector.
Definition: htransform.c:782
unsigned int GetAtomicNumber(unsigned int na, TMolecule *ml)
Returns the atomic number of a given atom.
Definition: babel.cpp:109
unsigned int nl
Definition: bioworld.h:37
void PrintWorld(char *fname, int argc, char **arg, Tworld *w)
Prints the world.
Definition: world.c:3902
void GetLinkTransformsFromSolutionPoint(Tparameters *p, boolean simp, double *sol, THTransform **tl, TLinkConf **def, Tworld *w)
Define transforms for the links from the a solution point.
Definition: world.c:1215
void DeleteID(void *a)
Destructor for identifiers.
Definition: vector.c:29
Relation between two links.
Definition: joint.h:168
void PrintBioWorld(Tparameters *p, char *fname, int argc, char **arg, TBioWorld *bw)
Prints the BioWorld information into a file.
Definition: bioworld.c:1919
void SetEquationType(unsigned int type, Tequation *eq)
Changes the type of the equation (SYSTEM_EQ, CARTESIAN_EQ, DUMMY_EQ, DERIVED_EQ). ...
Definition: equation.c:1076
#define NEW(_var, _n, _type)
Allocates memory space.
Definition: defines.h:385
void GenerateNormEquation(unsigned int vx, unsigned int vy, unsigned int vz, double n, Tequation *eq)
Construtor. Generates an equation that is the norm of a 3d vector.
Definition: equation.c:1494
Data structure to hold the information about the name of a file.
Definition: filename.h:271
unsigned int na
Definition: bioworld.h:31
void DeleteEquation(Tequation *eq)
Destructor.
Definition: equation.c:1859
void * GetVectorElement(unsigned int i, Tvector *vector)
Returns a pointer to a vector element.
Definition: vector.c:270
A homgeneous transform in R^3.
boolean IsAtomInProline(unsigned int na, TMolecule *ml)
Identifies atoms in proline residues.
Definition: babel.cpp:155
#define SYSTEM_VAR
One of the possible type of variables.
Definition: variable.h:24
unsigned int AddVariable2CS(Tvariable *v, TCuikSystem *cs)
Adds a variable to the system.
Definition: cuiksystem.c:2532
void BioWordSetAtomPositionsFromConformation(Tparameters *p, boolean simp, double *conformation, TBioWorld *bw)
Changes the position of the atoms.
Definition: bioworld.c:1852
void Atoms2Transforms(Tparameters *p, double *atoms, THTransform *t, TBioWorld *bw)
Generates a transform from a gobal frame to a local frame for each link.
Definition: bioworld.c:962
void GetAtomCoordinates(double *pos, TMolecule *ml)
Gets the positions of the atoms in the molecule.
Definition: babel.cpp:308
#define EQU
In a Tequation, the equation relational operator is equal.
Definition: equation.h:202
#define DUMMY_EQ
One of the possible type of equations.
Definition: equation.h:164
void DeleteCuikSystem(TCuikSystem *cs)
Destructor.
Definition: cuiksystem.c:5450
unsigned int GetCSSystemVars(boolean **sv, TCuikSystem *cs)
Identifies the system variables.
Definition: cuiksystem.c:2659
#define NORMAL_SHAPE
One of the possible type of shapes.
Definition: polyhedron.h:28
void WriteMolecule(char *fname, TMolecule *ml)
Writes the molecule information to a file.
Definition: babel.cpp:79
void ReadResidueList(char *fname, unsigned int *nr, char *ch, unsigned int **r)
Read the list of residues to be considered flexible.
Definition: bioworld.c:231
Definition of the Tfilename type and the associated functions.
void SetVariableInterval(Tinterval *i, Tvariable *v)
Sets the new range for the variable.
Definition: variable.c:68
#define TRUE
TRUE.
Definition: boolean.h:21
CBLAS_INLINE void Normalize(unsigned int s, double *v)
Normalizes a vector.
void NewSphere(double r, double *center, Tcolor *c, unsigned int gr, unsigned int bs, Tpolyhedron *p)
Constructor.
Definition: polyhedron.c:1020
void InitEquation(Tequation *eq)
Constructor.
Definition: equation.c:86
void Error(const char *s)
General error function.
Definition: error.c:80
unsigned int GetWorldNumSystemVariables(Tworld *w)
Number of system variables in the kinematic subsystem.
Definition: world.c:3297
void AddBond(unsigned int na1, unsigned int na2, TMolecule *ml)
Adds a bond between two atoms.
Definition: babel.cpp:209
#define NFUN
No trigonometric function for the variable.
Definition: variable_set.h:36
All the necessary information to generate equations for mechanisms.
Definition: world.h:229
A color.
Definition: color.h:86
void ReadRigidsAndHinges(char *fname, unsigned int **r, unsigned int *nh, unsigned int **h1, unsigned int **h2, boolean **added, TBioWorld *bw)
Read the list of rigids and hinges of a molecule.
Definition: bioworld.c:277
void DeleteLinkTransforms(THTransform *tl, TLinkConf *def, Tworld *w)
Deletes transforms for each link.
Definition: world.c:1336
void GetAtomPositions(char *fname, double *pos, TBioWorld *bw)
Initializes the atom positions in a BioWorld.
Definition: bioworld.c:437
Tvector joint1
Definition: bioworld.h:42
void SetEquationValue(double v, Tequation *eq)
Changes the right-hand value of the equation.
Definition: equation.c:1089
void AddMonomial(Tmonomial *f, Tequation *eq)
Adds a new monomial to the equation.
Definition: equation.c:1419
void TBondIterator
Iterator over the neighbours of a given atom.
Definition: babel.h:26
void DeleteWorld(Tworld *w)
Destructor.
Definition: world.c:3952
Tworld w
Definition: bioworld.h:29
A polyhedron.
Definition: polyhedron.h:134
void DeleteFileName(Tfilename *fn)
Destructor.
Definition: filename.c:205
unsigned int GetAtomResidue(unsigned int na, TMolecule *ml)
Returns the residue of a given atom.
Definition: babel.cpp:121
Tvector links
Definition: bioworld.h:38
void CheckAllCollisions(unsigned int fl, unsigned int fo, Tworld *w)
Activates all the possible collision between links and links and obstacles.
Definition: world.c:2097
unsigned int GetSolutionPointFromLinkTransforms(Tparameters *p, THTransform *tl, TLinkConf *def, double **sol, Tworld *w)
Determines the mechanisms configuration from the pose of all links.
Definition: world.c:1261
void HTransformFromVectors(double *x, double *y, double *z, double *h, THTransform *t)
Defines a homogeneous transform from 4 vectors.
Definition: htransform.c:341
void SetEquationCmp(unsigned int cmp, Tequation *eq)
Changes the relational operator (LEQ, GEQ, EQU) of the equation.
Definition: equation.c:1081
double Distance(unsigned int s, double *v1, double *v2)
Computes the distance of two points.
unsigned int nb
Definition: bioworld.h:32
void ResetMonomial(Tmonomial *f)
Reset the monomial information.
Definition: monomial.c:24
double DotProduct(double *v1, double *v2)
Computes the dot product of two 3d vectors.
Definition: geom.c:652
void HTransformInverse(THTransform *t, THTransform *ti)
Inverse of a homogeneous transform.
Definition: htransform.c:503
void RemoveBond(unsigned int na1, unsigned int na2, TMolecule *ml)
Removes a bond between two atoms.
Definition: babel.cpp:218
void DeleteBondIterator(TBondIterator *it)
Deletes the bond iterator defined in GetFirstNeighbour.
Definition: babel.cpp:299
double * localPos
Definition: bioworld.h:34
Tbox cutB
Definition: bioworld.h:52
void GenerateWorldEquations(Tparameters *p, Tworld *w)
Generates all the cuiksystems derived from the world information.
Definition: world.c:2431
An equation.
Definition: equation.h:237
void DeleteVariable(Tvariable *v)
Destructor.
Definition: variable.c:93
unsigned int GetNextNeighbour(unsigned int na, boolean *fix, TBondIterator *it, TMolecule *ml)
Gets the identifier of the next neighbour for a given atom in a molecule.
Definition: babel.cpp:271
void AdjustBioWorldGeometry(Tparameters *p, TBioWorld *bw)
Enforces all bond lengths and angles to be the same.
Definition: bioworld.c:1475
unsigned int AddLink2World(Tlink *l, boolean object, Tworld *w)
Adds a link to the mechanism in the world.
Definition: world.c:1383
void InitVector(unsigned int ele_size, void(*Copy)(void *, void *), void(*Delete)(void *), unsigned int max_ele, Tvector *vector)
Constructor.
Definition: vector.c:100
unsigned int BioWorldNAtoms(TBioWorld *bw)
Number of atoms in the molecule.
Definition: bioworld.c:1812
Definition of the Tjoint type and the associated functions.
void DifferenceVector(unsigned int s, double *v1, double *v2, double *v)
Substracts two vectors.
void NewRevoluteJoint(unsigned int id, unsigned int r, unsigned int linkID1, Tlink *link1, unsigned int linkID2, Tlink *link2, double **points, boolean hasLimits, Tinterval *range, double **rPoints, boolean avoidLimits, double avoidLimitsWeight, Tjoint *coupled, Tjoint *j)
Constructor.
Definition: joint.c:124
Minimalistic Cuik-OpenBabel interface.
void DetectLinksAndJointsFromResidues(unsigned int nr, char ch, unsigned int *rID, TBioWorld *bw)
Determines the rigid groups of atoms and the connections between them.
Definition: bioworld.c:553
void NewVariable(unsigned int type, char *name, Tvariable *v)
Constructor.
Definition: variable.c:20
char GetAtomChain(unsigned int na, TMolecule *ml)
Returns a char identifying the chain of the atom.
Definition: babel.cpp:138
TMolecule * ReadMolecule(char *fname)
Reads the molecule information from a file.
Definition: babel.cpp:46
boolean InitBioWorld(Tparameters *p, char *filename, unsigned int maxAtomsLink, double **conformation, TBioWorld *bw)
Initializes a world form a biomolecule.
Definition: bioworld.c:1394
void InitCuikSystem(TCuikSystem *cs)
Constructor.
Definition: cuiksystem.c:2167
void AddEquation2CS(Tparameters *p, Tequation *eq, TCuikSystem *cs)
Adds an equation to the system.
Definition: cuiksystem.c:2502
A scaled product of powers of variables.
Definition: monomial.h:32
A table of parameters.
void CreateFileName(char *path, char *name, char *suffix, char *ext, Tfilename *fn)
Constructor.
Definition: filename.c:22
void NewCylinder(double r, double *p1, double *p2, Tcolor *c, unsigned int gr, unsigned int bs, Tpolyhedron *p)
Constructor.
Definition: polyhedron.c:1044
boolean cut
Definition: bioworld.h:50
unsigned int AddJoint2World(Tjoint *j, Tworld *w)
Adds a joint to the mechanism in the world.
Definition: world.c:1436
A bridge between world structures and biological information.
double GetBoxDiagonal(boolean *used, Tbox *b)
Computes the diagonal of a (sub-)box.
Definition: box.c:678
Definition of a local chart on a manifold.
double * pos
Definition: bioworld.h:35
char * GetFileFullName(Tfilename *fn)
Gets the file full name (paht+name+extension).
Definition: filename.c:151
#define CT_REPRESENTATION
Representation.
Definition: parameters.h:215
void InitWorldFromMolecule(Tparameters *p, double **conformation, unsigned int maxAtomsLink, TBioWorld *bw)
Defines a world from molecular information.
Definition: bioworld.c:1130
A box.
Definition: box.h:83
Data associated with each variable in the problem.
Definition: variable.h:84
double VdWRadius(unsigned int na, TMolecule *ml)
Returns the Van der Waals radius for a given atom.
Definition: babel.cpp:172
Definition of the Tpolyhedron type and the associated functions.
#define MEM_DUP(_var, _n, _type)
Duplicates a previously allocated memory space.
Definition: defines.h:414
A cuiksystem, i.e., a set of variables and equations defining a position analysis problem...
Definition: cuiksystem.h:181
#define NO_UINT
Used to denote an identifier that has not been initialized.
Definition: defines.h:435
CBLAS_INLINE void SumVectorScale(unsigned int s, double *v1, double w, double *v2, double *v)
Adds two vectors with a scale.
Definition: basic_algebra.c:86
void InitBoxFromPoint(unsigned int dim, double *p, Tbox *b)
Initializes a box from a point.
Definition: box.c:43
THTransform * g2l
Definition: bioworld.h:47
void AddVariable2Monomial(unsigned int fn, unsigned int varid, unsigned int p, Tmonomial *f)
Adds a power variable to the monomial.
Definition: monomial.c:171
Tvector joint2
Definition: bioworld.h:43
unsigned int GetFirstNeighbour(unsigned int na, boolean *fix, TBondIterator **it, TMolecule *ml)
Gets the identifier of the first neighbour for a given atom in a molecule.
Definition: babel.cpp:244
void PrintCuikSystem(Tparameters *p, FILE *f, TCuikSystem *cs)
Prints a cuiksystem.
Definition: cuiksystem.c:5359
unsigned int BioWordConformationFromAtomPositions(Tparameters *p, double *atoms, double **conformation, TBioWorld *bw)
Produces the internal coordinates from the atom positions.
Definition: bioworld.c:1882
void GetCSVariableNames(char **varNames, TCuikSystem *cs)
Gets points to the variable names.
Definition: cuiksystem.c:2555
void DeleteBox(void *b)
Destructor.
Definition: box.c:1283
void PrintBoxSubset(FILE *f, boolean *used, char **varNames, Tbox *b)
Prints a (sub-)box.
Definition: box.c:1162
Tlink * GetWorldLink(unsigned int linkID, Tworld *w)
Gets a link from its identifier.
Definition: world.c:1842
Structure with the molecular and the mechanical information.
Definition: bioworld.h:28
#define ATOM_EXT
File extension for atom coordinates files.
Definition: filename.h:89
void HTransformDelete(THTransform *t)
Destructor.
Definition: htransform.c:887
unsigned int nj
Definition: bioworld.h:41
#define MEM_EXPAND(_var, _n, _type)
Expands a previously allocated memory space.
Definition: defines.h:404
double GetParameter(unsigned int n, Tparameters *p)
Gets the value for a particular parameter.
Definition: parameters.c:93
void InitWorld(Tworld *w)
Constructor.
Definition: world.c:1358
void DeletePolyhedron(Tpolyhedron *p)
Destructor.
Definition: polyhedron.c:1692
void ExpandBox(double *p, Tbox *b)
Expands a box so that it includes a given point.
Definition: box.c:67
void DeleteColor(Tcolor *c)
Destructor.
Definition: color.c:143
void NewInterval(double lower, double upper, Tinterval *i)
Constructor.
Definition: interval.c:47
void AddCt2Monomial(double k, Tmonomial *f)
Scales a monomial.
Definition: monomial.c:158
void SaveBioWorldBioInfo(Tparameters *p, char *fname, boolean simp, double *conformation, TBioWorld *bw)
Stores the BioWorld information in a molecular format (eg. pdb).
Definition: bioworld.c:1911
double BioWorldRMSE(double *pos, TBioWorld *bw)
Computes the RMSE.
Definition: bioworld.c:1860
#define HINGE_EXT
File extension for files storing the hinges of a molecule.
Definition: filename.h:245
Auxiliary functions to deal with sets of samples.
#define CT_VDW_RATIO
Ratio over over the Van der Waals radius.
Definition: parameters.h:408
unsigned int GetCSVariableID(char *name, TCuikSystem *cs)
Gets the numerical identifier of a variable given its name.
Definition: cuiksystem.c:2607
#define RIGID_EXT
File extension for files storing rigid part of molecules.
Definition: filename.h:237
void NewColor(double r, double g, double b, Tcolor *c)
Constructor.
Definition: color.c:14
void HTransformPrint(FILE *f, THTransform *t)
Prints the a homogeneous transform to a file.
Definition: htransform.c:822
double ComputeEnergy(TMolecule *ml)
Computes the potential energy of a molecule.
Definition: babel.cpp:334
Defines a interval.
Definition: interval.h:33
#define RES_EXT
File extension for files storing residue indentifiers.
Definition: filename.h:229
Tworld * BioWorldWorld(TBioWorld *bw)
Returns a pointer to the world generated from the bio-information.
Definition: bioworld.c:1807
void InitMonomial(Tmonomial *f)
Constructor.
Definition: monomial.c:17
void CrossProduct(double *v1, double *v2, double *v3)
Computes the cross product of two 3d vectors.
Definition: geom.c:639
void DetectLinksAndJointsFromRigidsAndHinges(unsigned int *rg, unsigned int nh, unsigned int *h1, unsigned int *h2, TBioWorld *bw)
Determines the rigid groups of atoms and the connections between them.
Definition: bioworld.c:766
void GetMoleculeBasicInfo(TBioWorld *bw)
Collects the basic information about the molecule.
Definition: bioworld.c:398
unsigned int * linkID
Definition: bioworld.h:39
void DeleteJoint(Tjoint *j)
Destructor.
Definition: joint.c:3190
#define DUMMY_VAR
One of the possible type of variables.
Definition: variable.h:53
unsigned int * nbAtom
Definition: bioworld.h:36
void DeleteMonomial(Tmonomial *f)
Destructor.
Definition: monomial.c:289
void GenerateDotProductEquation(unsigned int v1x, unsigned int v1y, unsigned int v1z, unsigned int v2x, unsigned int v2y, unsigned int v2z, unsigned int vc, double c, Tequation *eq)
Construtor. Generates the equation of the dot product of two unitary vectors.
Definition: equation.c:1588
void BioWordGetAtomPositionsFromConformation(Tparameters *p, boolean simp, double *conformation, double *pos, TBioWorld *bw)
Computes the position of the atoms.
Definition: bioworld.c:1822
unsigned int NewVectorElement(void *e, Tvector *vector)
Adds an element to the vector.
Definition: vector.c:216
#define DECOR_SHAPE
One of the possible type of shapes.
Definition: polyhedron.h:46
Link configuration.
Definition: link.h:276
void CopyID(void *a, void *b)
Copy constructor for identifiers.
Definition: vector.c:24
TMolecule * m
Definition: bioworld.h:30
unsigned int BioWorldConformationSize(TBioWorld *bw)
Number of variables used to represent a conformation.
Definition: bioworld.c:1817
#define HIDDEN_SHAPE
One of the possible type of shapes.
Definition: polyhedron.h:37