home *** CD-ROM | disk | FTP | other *** search
/ Nebula 1 / Nebula One.iso / Educational / MolViewer / Source / MolObj.m < prev    next >
Encoding:
Text File  |  1995-06-12  |  24.3 KB  |  1,037 lines
  1. /* MolObj.m - Copyright 1993  Steve Ludtke */
  2. /* This is the central object in the program. Almost everything */
  3. /* goes through this object on some level */
  4.  
  5. #import "Mtypes.h"
  6. #import "MolObj.h"
  7. #import "Inspector.h"
  8. #import "D3View.h"
  9. #import "AminoView.h"
  10.  
  11. #import <stdio.h>
  12. #import <string.h>
  13. #import <libc.h>
  14. #import <math.h>
  15. #import <appkit/appkit.h>
  16. #import <dpsclient/psops.h>
  17.  
  18. #define NSTEP 50
  19. #define NSTEP0 200
  20. #define swapa(n1,n2) { tmp=mol[0].am[n1]; mol[0].am[n1]=mol[0].am[n2]; mol[0].am[n2]=tmp; }
  21.  
  22. /* vector operations to simplify math, boy, c++ would be nice ... */
  23. float dot(float *a,float *b);
  24. void cross(float *a,float *b,float *c);
  25. void sub(float *a,float *b,float *c);
  26. float len(float *a);
  27. void count (Molecule *mol,int *cnt, int n1,short lev,short maxlev);
  28.  
  29. void tagatoms(Atom *atom,int n);
  30. float frand(float lo, float hi);
  31. extern ACID acid[AADEF];
  32.  
  33. @implementation MolObj
  34.  
  35. -init
  36. {
  37. FILE *in;
  38. char *bpath;
  39. char s[80];
  40. int i;
  41.  
  42. [super init];
  43.  
  44. srandom(time(0));
  45.  
  46. /* find where the app is located */
  47. bpath=(char *)[[NXBundle mainBundle] directory];
  48. /* read the covalent.rad file with atom preferences */
  49. sprintf(s,"%s/Library/MolViewer.dat",getenv("HOME"));
  50. in=fopen(s,"r");
  51. if (in==NULL) {
  52.     sprintf(s,"%s/MolViewer.dat",bpath);
  53.     in=fopen(s,"r");
  54.     if (in==NULL) { printf("No .rad\n"); exit(0); }
  55. }
  56. for (i=0; fgets(s,79,in)!=NULL&&i<NEL; i++) {
  57.     elinfo[i].color[0]=elinfo[i].color[1]=elinfo[i].color[2]=.7;
  58.     elinfo[i].name[0]=s[0];
  59.     if (s[1]!='\t') elinfo[i].name[1]=s[1];
  60.     else elinfo[i].name[1]=' ';
  61.     sscanf(&s[2],"%f %f %f %f",&elinfo[i].arad,&elinfo[i].color[0],&elinfo[i].
  62.         color[1],&elinfo[i].color[2]);
  63. /* the 1.3 is user adjustable after the program has been started */
  64.     elinfo[i].rad=1.3*elinfo[i].arad;
  65.     elinfo[i].image=nil;
  66. }
  67. fclose(in);
  68.  
  69. /* initialize arrays */
  70. for (i=0; i<MAXMOL; i++) { mol[i].verts=NULL; mol[i].nverts=NULL; }
  71.  
  72. /* load shapes for quick spacefilling mode */
  73. sprintf(s,"%s/h.tiff",bpath);
  74. elinfo[0].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
  75. sprintf(s,"%s/c.tiff",bpath);
  76. elinfo[5].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
  77. sprintf(s,"%s/n.tiff",bpath);
  78. elinfo[6].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
  79. sprintf(s,"%s/o.tiff",bpath);
  80. elinfo[7].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
  81. sprintf(s,"%s/p.tiff",bpath);
  82. elinfo[14].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
  83. sprintf(s,"%s/s.tiff",bpath);
  84. elinfo[15].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
  85.  
  86. nmol=0;
  87. stepmode=1;
  88. cstep=1.0;
  89. return self;
  90. }
  91.  
  92. /* we're the App's delegate, so tell everyone else that we're ready to go */
  93. -appDidInit:sender
  94. {
  95. [molView appDidInit:sender];
  96. [molView setData:mol :nmol :elinfo];
  97. [inspector startup:mol :nmol :elinfo];
  98. return self;
  99. }
  100.  
  101. /* read a MolViewer file (format in transition, don't use) */
  102. -readMol:sender
  103. {
  104. id panel;
  105. FILE *in;
  106. int i,j,k;
  107. float x,y,z,x1,y1,z1;
  108. char c[2],s[101];
  109.  
  110. x=y=z=0.0;
  111.  
  112. panel=[[OpenPanel new] allowMultipleFiles:NO];
  113. if (![panel runModalForTypes:NULL]) return self;
  114. in=fopen([panel filename],"r");
  115. if (in==NULL) return self;
  116.  
  117. [self remMol:self];
  118. while (fgets(s,100,in)!=NULL);
  119. if (sscanf(s," %d",&mol[nmol].numa)!=1 || mol[nmol].numa>5000)
  120.     { fclose(in); return self; }
  121. mol[nmol].numa++;
  122. mol[nmol].atom=malloc(sizeof(Atom)*(mol[nmol].numa+1));
  123. mol[nmol].coord=malloc(sizeof(RtPoint)*(mol[nmol].numa+1));
  124. rewind(in);
  125.  
  126. for (i=0; fgets(s,100,in)!=NULL; i++) {
  127.     if (s[0]=='#') { i=0; continue; }
  128.     sscanf(s,"%d %c%c %f %f %f",&j,c,&c[1],&x1,&y1,&z1);
  129.     if (c[1]=='\t') c[1]=' ';
  130.  
  131.     mol[nmol].coord[j][0]=x1;
  132.     mol[nmol].coord[j][1]=y1;
  133.     mol[nmol].coord[j][2]=z1;
  134.     mol[nmol].atom[j].numb=0;
  135.     mol[nmol].atom[j].cnum=i;    /* atom # within residue (for proteins) */
  136.     mol[nmol].atom[j].tag=mol[nmol].atom[j].sel=0;
  137.  
  138.     for (k=0; k<NEL; k++) if (c[0]==elinfo[k].name[0]&&c[1]==elinfo[k].name[1]) break;
  139.     if (k==NEL) printf("Can't find atom:%c\n",c[0]);
  140.     mol[nmol].atom[j].anum=k;
  141.  
  142.     x+=x1;    /* keep running total for finding average x,y,z */
  143.     y+=y1;
  144.     z+=z1;
  145. }
  146. fclose(in);
  147. j++;
  148. x/=(float)j;
  149. y/=(float)j;
  150. z/=(float)j;
  151.  
  152. /* move origin to center of molecule */
  153. for (i=0; i<mol[nmol].numa; i++) {
  154.     mol[nmol].coord[i][0]-=x;
  155.     mol[nmol].coord[i][1]-=y;
  156.     mol[nmol].coord[i][2]-=z;
  157. }
  158.  
  159. mol[nmol].numlb=0;
  160. strcpy(s,[panel filename]);
  161. s[strlen(s)-3]='b';
  162. s[strlen(s)-2]='n';
  163. s[strlen(s)-1]='d';
  164. in=fopen(s,"r");
  165. if (in==NULL) { printf("Can't open bnd\n"); nmol++; 
  166.     [molView setData:mol :nmol :elinfo]; return self; }
  167.  
  168. while (fgets(s,100,in)!=NULL) {
  169.     if (s[0]=='L') mol[nmol].numlb++;
  170. }
  171. mol[nmol].lb=malloc(sizeof(struct LBOND)*mol[nmol].numlb);
  172. rewind(in);
  173.  
  174. i=j=0;
  175. while (fgets(s,100,in)!=NULL) {
  176.     if (s[0]=='L') {
  177.         sscanf(&s[1],"%d %d %f %f",&mol[nmol].lb[i].n1,&mol[nmol].lb[i].n2, 
  178.             &mol[nmol].lb[i].d,&mol[nmol].lb[i].w);
  179.         if (mol[nmol].lb[i].n1>mol[nmol].numa-1 || 
  180.             mol[nmol].lb[i].n2>mol[nmol].numa-1) i--;
  181.         i++;
  182.     }
  183. }
  184. [self asnBnd:nmol];
  185. nmol++;
  186. [inspector newMol];
  187. [self update:U_TOTAL];
  188. return self;
  189. }
  190.  
  191. /* read a Protein Data Bank molecule */
  192. -readPdb:sender
  193. {
  194. id panel;
  195. FILE *in;
  196. int i,j,k,nb,b[8];
  197. float x,y,z;
  198. char s[101];
  199. RtPoint *a1,*a2;
  200.  
  201. panel=[[OpenPanel new] allowMultipleFiles:NO];
  202. if (![panel runModalForTypes:NULL]) return self;
  203. in=fopen([panel filename],"r");
  204. if (in==NULL) return self;
  205.  
  206. /* get rid of the current molecule */
  207. [self remMol:self];
  208. i=nb=0;
  209. k=1;
  210. x=y=z=0;
  211. /* first pass, count atoms and bonds */
  212. while (fgets(s,100,in)!=NULL) {
  213. /* still uses CONECT records to generate bonds. Needs to be replaced */
  214.     if (strncmp(s,"CONECT",6)==0) {
  215.         s[70]=0;
  216.         j=sscanf(&s[7]," %d %d %d %d %d %d %d %d",&k,&k,&k,&k,&k,&k,&k,&k);
  217.         nb+=j-1;
  218.     }
  219.     if (k && (strncmp(s,"ATOM",4)==0 || strncmp(s,"HETATM",6)==0)) {
  220.         sscanf(&s[7]," %d",&i);
  221.     }
  222. /* only reads the first structure in the file, TER terminates reading */
  223.     if (strncmp(s,"TER",3)==0) k=0;
  224. }
  225. rewind(in);
  226. mol[nmol].numa=i;
  227. mol[nmol].atom=malloc(sizeof(Atom)*(i+1));
  228. mol[nmol].coord=malloc(sizeof(RtPoint)*(i+1));
  229. mol[nmol].numlb=nb;
  230. if (nb!=0) mol[nmol].lb=malloc(sizeof(struct LBOND)*(nb+1));
  231.  
  232. i=nb=0;
  233. k=0;
  234. /* ok, this time actually read the data */
  235. while (fgets(s,100,in)!=NULL) {
  236.     if (strncmp(s,"CONECT",6)==0) {
  237.         s[70]=0;
  238.         j=sscanf(&s[7]," %d %d %d %d %d %d %d %d",b,b+1,b+2,b+3,b+4, 
  239.             b+5,b+6,b+7)-1;
  240.         for (i=0; i<j; i++) {
  241.             a1=(RtPoint *)&mol[nmol].coord[b[i]-1][0];
  242.             a2=(RtPoint *)&mol[nmol].coord[b[i+1]-1][0];
  243.             if (b[i]>mol[nmol].numa || b[i+1]>mol[nmol].numa) continue;
  244.             mol[nmol].lb[nb].n1=b[i]-1;
  245.             mol[nmol].lb[nb].n2=b[i+1]-1;
  246.             mol[nmol].lb[nb].w=1.0;
  247.             mol[nmol].lb[nb].d=sqrt(SQR(a1[0]-a2[0])+SQR(a1[1]-a2[1])+ 
  248.                 SQR(a1[2]-a2[2]));
  249.             nb++;
  250.         }
  251.     continue;
  252.     }
  253.     if (k) continue;
  254.     if (strncmp(s,"TER",3)==0) k=1;
  255.     if (strncmp(s,"ATOM",4)!=0 && strncmp(s,"HETATM",6)!=0) continue;
  256.     sscanf(&s[6]," %d",&i);
  257.     i--;
  258. /* try to find this atom in the atomic table */
  259.     for (j=0; j<NEL; j++) 
  260.         if(s[12]==elinfo[j].name[0]&&s[13]==elinfo[j].name[1]) break;
  261.     if (j==NEL) {
  262.         for (j=0; j<NEL; j++) if(s[13]==elinfo[j].name[0]) break;
  263.     }
  264.     if (j==NEL) printf("Can't find atom:%c%c\n",s[12],s[13]);
  265.     mol[nmol].atom[i].anum=j;
  266.     mol[nmol].atom[i].numb=0;
  267.     mol[nmol].atom[i].tag=mol[nmol].atom[i].sel=0;
  268.     sscanf(&s[30]," %f %f %f",&mol[nmol].coord[i][0], 
  269.         &mol[nmol].coord[i][1],&mol[nmol].coord[i][2]);
  270. /* find the "average" position for centering */
  271.     x+=mol[nmol].coord[i][0];
  272.     y+=mol[nmol].coord[i][1];
  273.     z+=mol[nmol].coord[i][2];
  274. }
  275. mol[nmol].numlb=nb;
  276.  
  277. x/=(float)mol[nmol].numa;
  278. y/=(float)mol[nmol].numa;
  279. z/=(float)mol[nmol].numa;
  280.  
  281. /* move origin to center of molecule */
  282. for (i=0; i<mol[nmol].numa; i++) {
  283.     mol[nmol].coord[i][0]-=x;
  284.     mol[nmol].coord[i][1]-=y;
  285.     mol[nmol].coord[i][2]-=z;
  286. }
  287.  
  288. [self asnBnd:nmol];
  289. nmol++;
  290. [inspector newMol];
  291. [self update:U_TOTAL];
  292. return self;
  293. }
  294.  
  295. /* read an Alchemy format file */
  296. -readAlch:sender
  297. {
  298. id panel;
  299. FILE *in;
  300. int i,j,k,l;
  301. float x,y,z,x1,y1,z1;
  302. float *a1,*a2;
  303. char c[4],s[101];
  304.  
  305. x=y=z=0.0;
  306.  
  307. panel=[[OpenPanel new] allowMultipleFiles:NO];
  308. if (![panel runModalForTypes:NULL]) return self;
  309. in=fopen([panel filename],"r");
  310. if (in==NULL) return self;
  311.  
  312. [self remMol:self];
  313. /* Ah Ha! Alchemy has a header with an atom/bond count */
  314. fgets(s,100,in);
  315. sscanf(s," %d",&mol[nmol].numa);
  316. sscanf(&s[13]," %d",&mol[nmol].numlb);
  317. mol[nmol].atom=malloc(sizeof(Atom)*(mol[nmol].numa+1));
  318. mol[nmol].coord=malloc(sizeof(RtPoint)*(mol[nmol].numa+1));
  319.  
  320. for (i=0; i<mol[nmol].numa; i++) {
  321.     fgets(s,100,in);
  322.     if (s[0]=='#') { i=0; continue; }
  323.     sscanf(s," %d %c%c%c%c %f %f %f",&j,c,&c[1],&c[2],&c[3],&x1,&y1,&z1);
  324.     j--;
  325.  
  326.     mol[nmol].coord[j][0]=x1;
  327.     mol[nmol].coord[j][1]=y1;
  328.     mol[nmol].coord[j][2]=z1;
  329.     mol[nmol].atom[j].numb=0;
  330.     mol[nmol].atom[j].cnum=i;    /* atom # within residue (for proteins) */
  331.  
  332. /* try to identify the atoms. Free bonds become hydrogens */
  333.     if (c[0]=='~') c[0]='H';
  334.     for (k=0; k<NEL; k++) {
  335.         if ((c[1]>='0'&&c[1]<='9')||c[1]==' '||c[2]!=' ') {
  336.             if (c[0]==elinfo[k].name[0]&&elinfo[k].name[1]==' ') break;
  337.         }
  338.         else if (c[0]==elinfo[k].name[0]&&c[1]==elinfo[k].name[1]) break;
  339.     }
  340.     if (k==NEL) printf("Can't find atom:%c\n",c[0]);
  341.     mol[nmol].atom[j].anum=k;
  342.     mol[nmol].atom[j].type[0]=c[1];
  343.     mol[nmol].atom[j].type[1]=c[2];
  344.     mol[nmol].atom[j].type[2]=c[3];
  345.     mol[nmol].atom[j].tag=mol[nmol].atom[j].sel=0;
  346.     x+=x1;    /* keep running total for finding average x,y,z */
  347.     y+=y1;
  348.     z+=z1;
  349. }
  350. j++;
  351. x/=(float)j;
  352. y/=(float)j;
  353. z/=(float)j;
  354.  
  355. /* move origin to center of molecule */
  356. for (i=0; i<mol[nmol].numa; i++) {
  357.     mol[nmol].coord[i][0]-=x;
  358.     mol[nmol].coord[i][1]-=y;
  359.     mol[nmol].coord[i][2]-=z;
  360. }
  361.  
  362. mol[nmol].lb=malloc(sizeof(struct LBOND)*mol[nmol].numlb);
  363.  
  364. /* now read the bonds */
  365. j=0;
  366. for (i=0; i<mol[nmol].numlb; i++) {
  367.     fgets(s,100,in);
  368.     sscanf(s,"%d %d %d",&j,&k,&l);
  369.     mol[nmol].lb[i].n1=k-1;
  370.     mol[nmol].lb[i].n2=l-1;
  371.     a1=(float *)&mol[nmol].coord[k-1][0];
  372.     a2=(float *)&mol[nmol].coord[l-1][0];
  373.     mol[nmol].lb[i].w=1.0;
  374.     mol[nmol].lb[i].type=s[19];
  375.     x=sqrt(SQR(a1[0]-a2[0])+SQR(a1[1]-a2[1])+SQR(a1[2]-a2[2]));
  376.     if (x>2.2) x=1.8;
  377.     mol[nmol].lb[i].d=x;
  378. }
  379.  
  380. /*for (i=0; i<mol[0].numa; i++) {
  381.     mol[0].coord[i][1]*=.5;
  382.     mol[0].coord[i][0]=frand(-4.0,4.0);
  383.     mol[0].coord[i][2]=frand(-4.0,4.0);
  384. }*/
  385.  
  386. [self asnBnd:nmol];
  387. nmol++;
  388. fclose(in);
  389. [inspector newMol];
  390. [self update:U_TOTAL];
  391. return self;
  392. }
  393.  
  394. /* delete the current molecule from memory */
  395. -remMol:sender
  396. {
  397. if (nmol==0) return self;
  398. nmol--;
  399. if (mol[nmol].numa!=0) free(mol[nmol].atom);
  400. if (mol[nmol].numlb!=0) free(mol[nmol].lb);
  401. if (mol[nmol].verts!=NULL) { free(mol[nmol].verts); mol[nmol].verts=NULL; }
  402. if (mol[nmol].nverts!=NULL) { free(mol[nmol].nverts); mol[nmol].nverts=NULL; }
  403. mol[nmol].numa=mol[nmol].numlb=0;
  404. [self update:U_TOTAL];
  405. return self;
  406. }
  407.  
  408. /* something has changed, so make sure all the displays are correct */
  409. -update:(int)level
  410. {
  411. if (level>=U_BONDS) [self setupMol];
  412. [molView setData:mol :nmol :elinfo];
  413. [inspector update:nmol :level];
  414. return self;
  415. }
  416.  
  417. /* generates a data table to speed up drawing */
  418. -setupMol
  419. {
  420. int i,j,k;
  421.  
  422. for (i=0; i<nmol; i++) {
  423.     k=mol[i].numlb;
  424.     if (mol[i].verts!=NULL) free(mol[i].verts);
  425.     if (mol[i].nverts!=NULL) free(mol[i].nverts);
  426.     mol[i].verts=malloc(sizeof(RtInt)*2*(k+1));
  427.     mol[i].nverts=malloc(sizeof(RtInt)*(k+1));
  428.     for (j=0; j<k; j++) {
  429.         mol[i].nverts[j]=2;
  430.         mol[i].verts[j*2]=mol[i].lb[j].n1;
  431.         mol[i].verts[j*2+1]=mol[i].lb[j].n2;
  432.     }
  433. }
  434. return self;
  435. }
  436.  
  437. /* Write an Alchemy format file */
  438. -dump:sender
  439. {
  440. int i;
  441. static id savePanel=nil;
  442. FILE *out;
  443. char s[80];
  444.  
  445. if (savePanel==nil) savePanel=[SavePanel new];
  446. [savePanel setRequiredFileType:"mol"];
  447.   
  448. if (![savePanel runModal]) return self;
  449. out=fopen([savePanel filename],"w");
  450.  
  451. fprintf(out,"%5d ATOMS, %5d BONDS,     0 CHARGES, UNK\n", 
  452.             mol[0].numa,mol[0].numlb);
  453. for (i=0; i<mol[0].numa; i++) {
  454.     fprintf(out,"%5d %c%c%c%c  %8.4f %8.4f %8.4f     0.0000\n", i+1,elinfo[mol[0].atom[i].anum].name[0],mol[0].atom[i].type[0],mol[0].atom[i].type[1],mol[0].atom[i].type[2],mol[0].coord[i][0],mol[0].coord[i][1], mol[0].coord[i][2]);
  455. }
  456.  
  457. for (i=0; i<mol[0].numlb; i++) {
  458.     if (mol[0].lb[i].type=='S') strcpy(s,"SINGLE");
  459.     else strcpy(s,"DOUBLE");
  460.     fprintf(out,"%5d %5d %5d  %s\n",i+1,mol[0].lb[i].n1+1,mol[0].lb[i].n2+1,s);
  461. }
  462. fclose(out);
  463. return self;
  464. }
  465.  
  466. /* copy bond table to redundant bond information in Atom structure */
  467. -asnBnd:(int)nm
  468. {
  469. int i,j;
  470. int lb,n1,n2;
  471. float l;
  472.  
  473. lb=mol[nm].numlb;
  474.  
  475. for (i=0; i<lb; i++) {
  476.     n1=mol[nm].lb[i].n1;
  477.     n2=mol[nm].lb[i].n2;
  478.     l=mol[nm].lb[i].d;
  479.  
  480.     j=mol[nm].atom[n1].numb++;
  481.     if (j>3) { j=3; /*printf("Too many bonds\n");*/ }
  482.     mol[nm].atom[n1].bnd[j]=n2;
  483.     mol[nm].atom[n1].bndl[j]=l;
  484.  
  485.     j=mol[nm].atom[n2].numb++;
  486.     if (j>3) { j=3; /*printf("Too many bonds\n");*/ }
  487.     mol[nm].atom[n2].bnd[j]=n1;
  488.     mol[nm].atom[n2].bndl[j]=l;
  489. }
  490.  
  491. for (i=0; i<mol[nm].numa; i++) {
  492.     for (j=mol[nm].atom[i].numb; j<4; j++) mol[nm].atom[i].bnd[j]=-1;
  493. }
  494.  
  495. if (nm==0) chi=[self calcChi:self];
  496. return self;
  497. }
  498.  
  499. /* locate and identify amino acids in the current molecule */
  500. -findAA:sender
  501. {
  502. int i,j,k,l,h=0,n,na;
  503. int cnt[16];
  504. struct AMINO tmp;
  505.  
  506. /* This is difficult to follow, but you'll just have to wade through it */
  507. /* if you need to change something. */
  508. na=0;
  509. for (i=0; i<mol[0].numa; i++) {
  510.     if (mol[0].atom[i].anum!=6 || mol[0].atom[i].type[0]!='A') continue;
  511.     mol[0].am[na].n=i;            /* assign N */
  512.     mol[0].atom[i].tag=1;
  513.     mol[0].am[na].cn= -1;
  514.                                 /* assign h,c1, and cn */
  515.     for (j=0; j<3; j++) {
  516.         k=mol[0].atom[i].bnd[j];
  517. /* h is attached to n */    
  518.         if (mol[0].atom[k].anum==0) mol[0].am[na].h=k;
  519. /* cn has an O attached, c1 has an H and the sidechain */
  520.         if (mol[0].atom[k].anum==5) {
  521.             if (mol[0].atom[k].type[0]=='3') mol[0].am[na].c1=k;
  522.             else if (mol[0].atom[k].type[0]=='2') mol[0].am[na].cn=k;
  523.             else printf("Bad carbon on Nam\n");
  524.         }
  525.     }
  526.                                 /* assign cn2 */
  527.     k=mol[0].am[na].cn;
  528.     if (k==-1) mol[0].am[na].cn2= -1;
  529.     else {
  530. /* cn2 is attached to cn and an o */
  531.         for (j=0; j<mol[0].atom[k].numb; j++) {
  532.             n=mol[0].atom[k].bnd[j];
  533.             if (mol[0].atom[n].anum==5 && mol[0].atom[n].type[0]=='3') 
  534.                 mol[0].am[na].cn2=n;
  535.         }
  536.     }
  537.                                 /* assign c2 */
  538.     k=mol[0].am[na].c1;
  539.     mol[0].atom[k].tag=1;
  540. /* c2 is attached to c1, h, and a sidechain */
  541.     for (j=0; j<mol[0].atom[k].numb; j++) {
  542.         n=mol[0].atom[k].bnd[j];
  543.         if (mol[0].atom[n].anum==5 && mol[0].atom[n].type[0]=='2') 
  544.             mol[0].am[na].c2=n;
  545.     }
  546.                                 /* assign o */
  547.     k=mol[0].am[na].c2;
  548.     for (j=0; j<mol[0].atom[k].numb; j++) {
  549.         n=mol[0].atom[k].bnd[j];
  550.         if (mol[0].atom[n].anum==7) mol[0].am[na].o=n;
  551.     }
  552.                                 /* assign sc */
  553.     k=mol[0].am[na].c1;
  554.     for (j=0; j<4; j++) {
  555.         n=mol[0].atom[k].bnd[j];
  556.         if (n==mol[0].am[na].c2||n==mol[0].am[na].cn) continue;
  557.         if (mol[0].atom[n].anum==0) { h=n; continue; }
  558.         if (mol[0].atom[n].anum==5) { mol[0].am[na].sc=n; break; }
  559.     }
  560. /* if glycine, an h is chosen randomly */    
  561.     if (j==4) mol[0].am[na].sc=h;
  562.  
  563.     na++;
  564.     if (na==MAXAA) { na--; printf("Max # of amino acids exceeded\n"); }
  565. }    
  566.  
  567. mol[0].numaa=na;
  568.  
  569. for (i=0; i<mol[0].numaa; i++) {
  570. /* calculate the dihedrals */
  571.     mol[0].am[i].sel=0;
  572.     mol[0].am[i].omega=[self dihedral:mol[0].am[i].c1 :mol[0].am[i].n :mol[0].am[i].cn :mol[0].am[i].cn2];
  573.     mol[0].am[i].psi=[self dihedral:mol[0].am[i].n :mol[0].am[i].c1 :mol[0].am[i].c2 :mol[0].am[i].o]-M_PI;
  574.     mol[0].am[i].phi=[self dihedral:mol[0].am[i].c2 :mol[0].am[i].c1 :mol[0].am[i].n :mol[0].am[i].h]-M_PI;
  575.     mol[0].am[i].anum=0;
  576.     if (mol[0].am[i].psi< -M_PI) mol[0].am[i].psi+=M_PI;
  577.     if (mol[0].am[i].phi< -M_PI) mol[0].am[i].phi+=M_PI;
  578.  
  579. /* try to identify the individual amino acids by counting and a little */
  580. /* chain testing when the result is ambigous. Exotic amino acids */
  581. /* might be incorrectly assigned. Proline is also incorrectly assigned. */ 
  582.     for (j=0; j<16; j++) cnt[j]=0;
  583.     j=mol[0].am[i].sc;
  584.     if (mol[0].atom[j].anum==0) mol[0].am[i].anum=9;
  585.     else {
  586.         count(mol,cnt,mol[0].am[i].sc,0,15);    
  587.         for (k=0; k<AADEF; k++) {
  588.             if (acid[k].nc==cnt[5]&&acid[k].ns==cnt[15]&&acid[k].no==cnt[7]&& acid[k].nn==cnt[6]) break;
  589.         }
  590.         if (k==AADEF) { mol[0].am[i].anum=0;
  591.             printf("C=%d S=%d O=%d N=%d\n",cnt[5],cnt[15],cnt[7],cnt[6]);
  592.         }
  593.         else mol[0].am[i].anum=k;
  594.     }
  595.     if     (mol[0].am[i].anum==3) {
  596.         l=mol[0].am[i].sc;
  597.         cnt[5]=0;
  598.         for (k=0; k<mol[0].atom[l].numb; k++) if (mol[0].atom[mol[0].atom[l].bnd[k]].anum==5) cnt[5]++;
  599.         if (cnt[5]==3) mol[0].am[i].anum=4;
  600.     }
  601. }
  602. for (i=0; i<mol[0].numa; i++) mol[0].atom[i].tag=0;
  603.  
  604. /* sort amino acids so N terminus is FIRST */
  605. for (i=1; i<mol[0].numaa; i++) if (mol[0].am[i].cn==-1) swapa(0,i);
  606. for (i=1; i<mol[0].numaa; i++) {
  607.     for (j=i+2; j<mol[0].numaa; j++) {
  608.         if (mol[0].am[i-1].c2==mol[0].am[j].cn) swapa(i,j);
  609.     }
  610. }
  611. return self;
  612. }
  613.  
  614. /* recursive routine to decend the sidechain and count the number of */
  615. /* atoms of each type */
  616. void count (Molecule *mol,int *cnt, int n1,short lev,short maxlev)
  617. {
  618. int i,j;
  619.  
  620. if (lev>maxlev||mol[0].atom[n1].tag) return;
  621. j=mol[0].atom[n1].anum;
  622. if (j<16) cnt[j]++;
  623. mol[0].atom[n1].tag=1;
  624. for (i=0; i<mol[0].atom[n1].numb; i++) {
  625.     count(mol,cnt,mol[0].atom[n1].bnd[i],lev+1,maxlev);
  626. }
  627. return;
  628. }
  629.  
  630. /* calculates a "chi" for fitting (not used yet ...) */
  631. -(float)calcChi:sender
  632. {
  633. int i;
  634. int lb;
  635. float l,ch;
  636. float *a1,*a2;
  637.  
  638. lb=mol[0].numlb;
  639. ch=0.0;
  640. for (i=0; i<lb; i++) {
  641.     a1=(float *)&mol[0].coord[mol[0].lb[i].n1];
  642.     a2=(float *)&mol[0].coord[mol[0].lb[i].n2];
  643.     l=mol[0].lb[i].d;
  644.     ch+=fabs(l-sqrt(SQR(a1[0]-a2[0])+SQR(a1[1]-a2[1])+SQR(a1[2]-a2[2])));
  645. }
  646. ch/=lb;
  647. return ch;
  648. }
  649.  
  650. /* fitting stuff ... */
  651. -minStep:sender
  652. {
  653. int i,j,k,m,c,c2;
  654. float x,y,z,d,d2,v1,v2;
  655. float x1,y1,z1;
  656. float *a0,*a1;
  657. short *ba;
  658.  
  659. c2=0;
  660. if (stepmode==0) {
  661.     for (i=0; i<NSTEP0; i++) {
  662.         x1=y1=z1=0.0;
  663.         j=random()%mol[0].numa;
  664.         k=mol[0].atom[j].numb;
  665.         a0=(float *)&mol[0].coord[j];
  666.         for (m=0; m<k; m++) {
  667.             a1=(float *)&mol[0].coord[mol[0].atom[j].bnd[m]];
  668.             x=a1[0]-a0[0];
  669.             y=a1[1]-a0[1];
  670.             z=a1[2]-a0[2];
  671.             d2=sqrt(SQR(x)+SQR(y)+SQR(z));
  672.             d=(1.0/mol[0].atom[j].bndl[m]-1.0/d2);
  673.             x1+=x*d;
  674.             y1+=y*d;
  675.             z1+=z*d;
  676.         }
  677.         x1/=(float)k;
  678.         y1/=(float)k;
  679.         z1/=(float)k;
  680.         a0[0]+=x1;
  681.         a0[1]+=y1;
  682.         a0[2]+=z1;
  683.     }
  684.     chi=[self calcChi:self];
  685.     c2=100;
  686. }
  687. else if (stepmode==1) {
  688.     for (i=0; i<NSTEP; i++) {
  689.         j=random()%mol[0].numa;
  690.         m=1;
  691.         c=0;
  692.  
  693.         ba=mol[0].atom[j].bnd;
  694.         a0=(float *)&mol[0].coord[j];
  695.         v1=0.0;
  696.         for (k=0; k<mol[0].numa; k++) {
  697.             if (k==j||k==ba[0]||k==ba[1]||k==ba[2]||k==ba[3]) continue;
  698.             a1=(float *)&mol[0].coord[k];
  699.             d=sqrt(SQR(a0[0]-a1[0])+SQR(a0[1]-a1[1])+SQR(a0[2]-a1[2]));
  700.             v1+=1.0/d;
  701.         }
  702.  
  703.         while (m&&c++<30) {
  704.             x1=frand(-cstep,cstep);
  705.             y1=frand(-cstep,cstep);
  706.             z1=frand(-cstep,cstep);
  707.  
  708.             v2=0.0;
  709.             a0[0]+=x1;
  710.             a0[1]+=y1;
  711.             a0[2]+=z1;
  712.  
  713.             m=0;            
  714.             for (k=0; k<mol[0].numa; k++) {
  715.                 if (k==j||k==ba[0]||k==ba[1]||k==ba[2]||k==ba[3]) continue;
  716.                 a1=(float *)&mol[0].coord[k];
  717.                 d=sqrt(SQR(a0[0]-a1[0])+SQR(a0[1]-a1[1])+SQR(a0[2]-a1[2]));
  718.                 v2+=1.0/d;
  719.             }
  720. /*            if (v1<v2) m=1; else m=0;
  721.             if (m) { a0[0]-=x1; a0[1]-=y1; a0[2]-=z1; continue; } */
  722.             d=v1/(mol[0].numa*6.0);
  723.             v2/=mol[0].numa*6.0;
  724.             x=[self calcChi:self];
  725.             if (x+v2>chi+d) { m=1; a0[0]-=x1; a0[1]-=y1; a0[2]-=z1; continue; }
  726.             chi=x;
  727.         }
  728.         if (c<=30) c2++;
  729.     }
  730.     [chiDevDis setFloatValue:chi];
  731. }
  732. [c2Dis setIntValue:c2];
  733. if (c2<3) {
  734.     cstep*=.9;
  735.     [csDis setFloatValue:cstep];
  736. }
  737. return self;
  738. }
  739.  
  740. /* more fitting stuff */
  741. -resetCs:sender
  742. {
  743. cstep=1.0;
  744. [csDis setFloatValue:cstep];
  745. return self;
  746. }
  747. /* this too ... */
  748. -togFit:sender
  749. {
  750. stepmode=[sender intValue];
  751. return self;
  752. }
  753.  
  754. /* test routine (not used) */
  755. -rotest:sender
  756. {
  757. printf("1: %f\n",[self dihedral:5 :1 :3 :6]*DRC);
  758. [self rotateAbout:.0873 :1 :3];
  759. printf("2: %f\n",[self dihedral:5 :1 :3 :6]*DRC);
  760. [self update:U_POS];
  761. return self;
  762. }
  763.  
  764. /* set amino acid dihedral angles */
  765. -setOmega:(float)ang
  766. {
  767. int i;
  768. float a,b;
  769.  
  770. a=ang*DRC;
  771. for (i=0; i<mol[0].numaa; i++) {
  772.     if (mol[0].am[i].cn==-1||(!mol[0].am[i].sel)) continue;
  773.     mol[0].am[i].omega=a;
  774.     b=[self dihedral:mol[0].am[i].c1 :mol[0].am[i].n :mol[0].am[i].cn
  775.         :mol[0].am[i].cn2];
  776.     [self rotateAbout:a-b :mol[0].am[i].n :mol[0].am[i].cn];    
  777. }
  778. [self centerMol:0];
  779. return self;
  780. }
  781.  
  782. -setPsi:(float)ang
  783. {
  784. int i;
  785. float a,b;
  786.  
  787. a=ang+180.0;
  788. if (a>180.0) a-=360.0;
  789. a*=DRC;
  790.  
  791. for (i=0; i<mol[0].numaa; i++) {
  792.     if (!mol[0].am[i].sel) continue;
  793.     mol[0].am[i].psi=ang*DRC;
  794.     b=[self dihedral:mol[0].am[i].n :mol[0].am[i].c1 :mol[0].am[i].c2
  795.         :mol[0].am[i].o];
  796.     [self rotateAbout:a-b :mol[0].am[i].c1 :mol[0].am[i].c2];    
  797. }
  798. [self centerMol:0];
  799. return self;
  800. }
  801.  
  802. -setPhi:(float)ang
  803. {
  804. int i;
  805. float a,b;
  806.  
  807. a=ang+180.0;
  808. if (a>180.0) a-=360.0;
  809. a*=DRC;
  810.  
  811. for (i=0; i<mol[0].numaa; i++) {
  812.     if (!mol[0].am[i].sel) continue;
  813.     mol[0].am[i].phi=ang*DRC;
  814.     b=[self dihedral:mol[0].am[i].c2 :mol[0].am[i].c1 :mol[0].am[i].n
  815.         :mol[0].am[i].h];
  816.     [self rotateAbout:a-b :mol[0].am[i].c1 :mol[0].am[i].n];    
  817. }
  818. [self centerMol:0];
  819. return self;
  820. }
  821.  
  822. /* take 4 atoms and calculate a dihedral angle. The math is right, but */
  823. /* it would take too much space to explain it. */
  824. -(float)dihedral:(int)n1 :(int)n2 :(int)n3 :(int)n4
  825. {
  826. float cb[3],ca[3],bd[3],bc[3],v1[3],v2[3],v3[3];
  827. float t,t2;
  828.  
  829. sub(&mol[0].coord[n2][0],&mol[0].coord[n3][0],cb);
  830. sub(&mol[0].coord[n3][0],&mol[0].coord[n2][0],bc);
  831. sub(&mol[0].coord[n1][0],&mol[0].coord[n3][0],ca);
  832. sub(&mol[0].coord[n4][0],&mol[0].coord[n2][0],bd);
  833.  
  834. cross(cb,ca,v1);
  835. cross(bd,bc,v2);
  836. cross(v1,v2,v3);
  837. t=dot(v3,bc)/len(bc);
  838. t2=dot(v1,v2);
  839.  
  840. return(atan2(t,t2));
  841. }
  842.  
  843. /* calculate a bond angle */
  844. -(float)bondAng:(int)n1 :(int)n2 :(int)n3;
  845. {
  846. float v1[3],v2[3],r;
  847.  
  848. sub(&mol[0].coord[n1][0],&mol[0].coord[n2][0],v1);
  849. sub(&mol[0].coord[n3][0],&mol[0].coord[n2][0],v2);
  850.  
  851. r=dot(v1,v2)/(len(v1)*len(v2));
  852. return (acos(r));
  853. }
  854.  
  855. /* center the molecule on the screen */
  856. - centerMol:(int)m
  857. {
  858. float x,y,z;
  859. int i,j;
  860.  
  861. x=y=z=0.0;
  862. for (i=0; i<mol[m].numa; i++) {
  863.     x+=mol[m].coord[i][0];
  864.     y+=mol[m].coord[i][1];
  865.     z+=mol[m].coord[i][2];
  866. }
  867.  
  868. x/=(float)mol[m].numa;
  869. y/=(float)mol[m].numa;
  870. z/=(float)mol[m].numa;
  871.  
  872. for (i=0; i<mol[m].numa; i++) {
  873.     mol[m].coord[i][0]-=x;
  874.     mol[m].coord[i][1]-=y;
  875.     mol[m].coord[i][2]-=z;
  876. }
  877.  
  878. i=mol[m].am[4].o;
  879. j=mol[m].am[7].h;
  880. x=sqrt(SQR(mol[m].coord[i][0]-mol[m].coord[j][0])+
  881.     SQR(mol[m].coord[i][1]-mol[m].coord[j][1])+
  882.     SQR(mol[m].coord[i][2]-mol[m].coord[j][2]));
  883.  
  884. i=mol[m].am[4].o;
  885. j=mol[m].am[6].h;
  886. y=sqrt(SQR(mol[m].coord[i][0]-mol[m].coord[j][0])+
  887.     SQR(mol[m].coord[i][1]-mol[m].coord[j][1])+
  888.     SQR(mol[m].coord[i][2]-mol[m].coord[j][2]));
  889.  
  890. [hba setFloatValue:x];
  891. [hb3 setFloatValue:y];
  892.  
  893. [self update:U_POS];
  894. return self;
  895. }
  896.  
  897. /* set the view window to a variety of convenient sizes */
  898. -set320:sender
  899. {
  900. [vWindow sizeWindow:324 :202];
  901. return self;
  902. }
  903.  
  904. -set640:sender
  905. {
  906. [vWindow sizeWindow:644 :402];
  907. return self;
  908. }
  909.  
  910. -set6404:sender
  911. {
  912. [vWindow sizeWindow:644 :482];
  913. return self;
  914. }
  915.  
  916. -set533:sender
  917. {
  918. [vWindow sizeWindow:537 :402];
  919. return self;
  920. }
  921.  
  922. /* do a dihedral rotation about the line connecting n1 and n2 */
  923. -rotateAbout:(float)a :(int)n1 :(int)n2
  924. {
  925. Atom *atom;
  926. int i;
  927. RtPoint n,r;
  928. float mx[9],m;
  929.  
  930. a= -a;
  931. atom=mol[0].atom;
  932. for (i=0; i<mol[0].numa; i++) atom[i].tag=0;
  933. atom[n1].tag=1;
  934. tagatoms(atom,n2);
  935. atom[n1].tag=0;
  936. atom[n2].tag=0;
  937.  
  938. n[0]=mol[0].coord[n2][0]-mol[0].coord[n1][0]; /* n==vec from n1 to n2 */
  939. n[1]=mol[0].coord[n2][1]-mol[0].coord[n1][1];
  940. n[2]=mol[0].coord[n2][2]-mol[0].coord[n1][2];
  941. m=sqrt(SQR(n[0])+SQR(n[1])+SQR(n[2]));
  942. n[0]/=m;    /* n is unit vector */
  943. n[1]/=m;
  944. n[2]/=m;
  945.  
  946. m=1.0-cos(a);    /* set up mx transformation matrix */
  947. mx[0]=cos(a)+n[0]*n[0]*m;
  948. mx[1]=n[0]*n[1]*m+n[2]*sin(a);
  949. mx[2]=n[0]*n[2]*m-n[1]*sin(a);
  950. mx[3]=n[0]*n[1]*m-n[2]*sin(a);
  951. mx[4]=cos(a)+n[1]*n[1]*m;
  952. mx[5]=n[2]*n[1]*m+n[0]*sin(a);
  953. mx[6]=n[0]*n[2]*m+n[1]*sin(a);
  954. mx[7]=n[1]*n[2]*m-n[0]*sin(a);
  955. mx[8]=cos(a)+n[2]*n[2]*m;
  956.  
  957. for (i=0; i<mol[0].numa; i++) {
  958.     if (atom[i].tag==0) continue;
  959.     r[0]=mol[0].coord[i][0]-mol[0].coord[n1][0];
  960.     r[1]=mol[0].coord[i][1]-mol[0].coord[n1][1];
  961.     r[2]=mol[0].coord[i][2]-mol[0].coord[n1][2];
  962.  
  963.     mol[0].coord[i][0]=r[0]*mx[0]+r[1]*mx[1]+r[2]*mx[2]+mol[0].coord[n1][0];
  964.     mol[0].coord[i][1]=r[0]*mx[3]+r[1]*mx[4]+r[2]*mx[5]+mol[0].coord[n1][1];
  965.     mol[0].coord[i][2]=r[0]*mx[6]+r[1]*mx[7]+r[2]*mx[8]+mol[0].coord[n1][2];
  966.     atom[i].tag=0;
  967. }
  968.  
  969. return self;
  970. }
  971.  
  972. /* select atoms from SelectView data Structure */
  973. -setSelAtom:(struct SELDAT *)list
  974. {
  975. int i;
  976.  
  977. for (i=0; i<mol[0].numa; i++) {
  978.     if (list[i].sel) mol[0].atom[i].sel=1;
  979.     else mol[0].atom[i].sel=0;
  980. }
  981.  
  982. [molView display];
  983. return self;
  984. }
  985.  
  986. /* recursive atom tagger for rotation routines */
  987. void tagatoms(Atom *atom,int n)
  988. {
  989. int i,j;
  990.  
  991. for (i=0; i<atom[n].numb; i++) {
  992.     j=atom[n].bnd[i];
  993.     if (atom[j].tag==0) {
  994.         atom[j].tag=1;
  995.         tagatoms(atom,j);
  996.     }
  997. }
  998. return;
  999. }
  1000.  
  1001. /* dot product */
  1002. float dot(float *a,float *b)
  1003. {
  1004. return (a[0]*b[0]+a[1]*b[1]+a[2]*b[2]);
  1005. }
  1006.  
  1007. /* cross product */
  1008. void cross(float *a,float *b,float *c)
  1009. {
  1010. c[0]=a[1]*b[2]-a[2]*b[1];
  1011. c[1]=a[2]*b[0]-a[0]*b[2];
  1012. c[2]=a[0]*b[1]-a[1]*b[0];
  1013. return;
  1014. }
  1015.  
  1016. /* vector subtract */
  1017. void sub(float *a,float *b,float *c)
  1018. {
  1019. c[0]=a[0]-b[0];
  1020. c[1]=a[1]-b[1];
  1021. c[2]=a[2]-b[2];
  1022. return;
  1023. }
  1024.  
  1025. /* vector length */
  1026. float len(float *a)
  1027. {
  1028. return (sqrt(SQR(a[0])+SQR(a[1])+SQR(a[2])));
  1029. }
  1030.  
  1031. /* return a random float between lo and hi */
  1032. float frand(float lo, float hi)
  1033. {
  1034.     return (((float)random() / (float)MAXINT) * (hi - lo) + lo);
  1035. }
  1036.  
  1037. @end