This is MolObj.m in view mode; [Download] [Up]
/* MolObj.m - Copyright 1993 Steve Ludtke */ /* This is the central object in the program. Almost everything */ /* goes through this object on some level */ #import "Mtypes.h" #import "MolObj.h" #import "Inspector.h" #import "D3View.h" #import "AminoView.h" #import "NrgObj.h" #import <stdio.h> #import <string.h> #import <libc.h> #import <math.h> #import <ctype.h> #import <appkit/appkit.h> #import <dpsclient/psops.h> #define NSTEP 50 #define NSTEP0 200 #define swapa(n1,n2) { tmp=mol[0].am[n1]; mol[0].am[n1]=mol[0].am[n2]; mol[0].am[n2]=tmp; } /* vector operations to simplify math, boy, c++ would be nice ... */ float dot(float *a,float *b); void cross(float *a,float *b,float *c); void sub(float *a,float *b,float *c); float len(float *a); void count (Molecule *mol,int *cnt, int n1,short lev,short maxlev); #define mul(a,b) { a[0]*=b; a[1]*=b; a[2]*=b; } #define add2(a,b) { a[0]+=b[0]; a[1]+=b[1]; a[2]+=b[2]; } #define sub2(a,b) { a[0]-=b[0]; a[1]-=b[1]; a[2]-=b[2]; } #define add(a,b,c) { c[0]=a[0]+b[0]; c[1]=a[1]+b[1]; c[2]=a[2]+b[2]; } void tagatoms(Atom *atom,int n); float frand(float lo, float hi); extern ACID acid[AADEF]; @implementation MolObj -init { char *bpath; char s[140]; int i; [super init]; srandom(time(0)); /* find where the app is located */ bpath=(char *)[[NXBundle mainBundle] directory]; /* initialize arrays */ for (i=0; i<MAXMOL; i++) { mol[i].verts=NULL; mol[i].nverts=NULL; } sprintf(s,"%s/Library/MolViewer",getenv("HOME")); if (access(s,0)!=0) mkdir(s,0xffff); /* load shapes for quick spacefilling mode */ for (i=0; i<NEL; i++) { elinfo[i].image=nil; sprintf(s,"%s/Library/MolViewer/el%02d.tiff",getenv("HOME"),i); if (access(s,0)) continue; elinfo[i].image=[[[NXImage alloc] initFromFile:s] setScalable:YES]; } elinfo[0].name[0]='0'; /* flag to see if appDidInit */ nmol=0; mol[0].numa=mol[0].numlb=0; stepmode=1; cstep=1.0; frag.atom=NULL; updFlag=1; return self; } /* we're the App's delegate, so tell everyone else that we're ready to go */ -appDidInit:sender { FILE *in; char *bpath; char s[140]; int i; /* find where the app is located */ bpath=(char *)[[NXBundle mainBundle] directory]; /* read the MolViewer.dat file with atom preferences */ sprintf(s,"%s/Library/MolViewer/MolViewer.dat",getenv("HOME")); in=fopen(s,"r"); if (in==NULL) { sprintf(s,"open %s/README.rtf",bpath); system(s); sprintf(s,"%s/MolViewer.dat",bpath); in=fopen(s,"r"); if (in==NULL) [self error:"No MolViewer.dat file!" :YES]; } for (i=0; fgets(s,79,in)!=NULL&&i<NEL; i++) { elinfo[i].color[0]=elinfo[i].color[1]=elinfo[i].color[2]=.7; elinfo[i].name[0]=s[0]; if (s[1]!='\t') elinfo[i].name[1]=s[1]; else elinfo[i].name[1]=' '; sscanf(&s[2],"%f %f %f %f",&elinfo[i].arad,&elinfo[i].color[0],&elinfo[i]. color[1],&elinfo[i].color[2]); /* the 1.3 is user adjustable after the program has been started */ elinfo[i].rad=1.3*elinfo[i].arad; } if (i<NEL) { /* old file, add 2 new elms */ [self error:"Old MolViewer.dat file found. Please save new one." :NO]; elinfo[84].name[0]='*'; elinfo[84].name[1]=' '; elinfo[84].color[0]=elinfo[84].color[1]=elinfo[84].color[2]=1.0; elinfo[84].arad=.2; elinfo[84].rad=.26; elinfo[85].name[0]='L'; elinfo[85].name[1]=' '; elinfo[85].color[0]=elinfo[84].color[1]=elinfo[84].color[2]=.4; elinfo[85].arad=.2; elinfo[85].rad=.26; } fclose(in); [molView appDidInit:sender]; [molView setData:mol :nmol :elinfo]; [inspector startup:mol :nmol :elinfo]; [nrgobj startup:mol :elinfo]; [self app:self openFile:NULL type:NULL]; return self; } -(BOOL)appAcceptsAnotherFile:sender { return YES; } /* open a file from the workspace */ -(int)app:sender openFile:(const char *)file type:(const char *)ftype { static char s[120] = { 0 },t[10]; if (elinfo[0].name[0]=='0') { strcpy(s,file); strcpy(t,ftype); return YES; } if (file==NULL) f2open=(void *)s; else { f2open=(void *)file; strcpy(t,ftype); } if (strlen(f2open)==0) return NO; if (strcmp(t,"mvb")==0||strcmp(t,"mvt")==0) [self readMV:nil]; else if (strcmp(t,"mol")==0) [self readAlch:nil]; else if (strcmp(t,"ent")==0) [self readPdb:nil]; else return NO; return YES; } /* read a MolViewer file */ -readMV:sender { id panel; FILE *in; int i; struct MVatm atm; struct MVbnd bnd; unsigned short n; if (sender==nil) in=fopen(f2open,"r"); else { panel=[[OpenPanel new] allowMultipleFiles:NO]; if (![panel runModalForTypes:NULL]) return self; in=fopen([panel filename],"r"); } if (in==NULL) return self; fread(&i,4,1,in); if (strncmp((char *)&i,"MVTX",4)==0) { [self readMVTX:in]; return self; } if (strncmp((char *)&i,"MVB",4)) { fclose(in); return self; } [self remMol:self]; fread(&n,2,1,in); mol[nmol].numa=n; fread(&n,2,1,in); mol[nmol].numlb=n; fread(&i,4,1,in); /* for future expansion */ mol[nmol].atom=malloc(sizeof(Atom)*(mol[nmol].numa+1)); mol[nmol].coord=malloc(sizeof(RtPoint)*(mol[nmol].numa+1)); mol[nmol].lb=malloc(sizeof(struct LBOND)*mol[nmol].numlb); for (i=0; i<mol[nmol].numa; i++) { fread(&atm,18,1,in); if (atm.anum==255) atm.anum=84; mol[nmol].atom[i].anum=atm.anum; if (atm.type<32) atm.type=32; mol[nmol].atom[i].type=atm.type; mol[nmol].coord[i][0]=atm.x; mol[nmol].coord[i][1]=atm.y; mol[nmol].coord[i][2]=atm.z; mol[nmol].atom[i].charge=atm.c; mol[nmol].atom[i].numb=0; mol[nmol].atom[i].tag=mol[nmol].atom[i].sel=mol[nmol].atom[i].lab=0; } for (i=0; i<mol[nmol].numlb; i++) { fread(&bnd,5,1,in); mol[nmol].lb[i].n1=bnd.n1; mol[nmol].lb[i].n2=bnd.n2; mol[nmol].lb[i].type=bnd.type; mol[nmol].lb[i].d=mol[nmol].lb[i].w=0; } mol[nmol].ell=NULL; fclose(in); [self asnBnd:nmol]; nmol++; [inspector newMol]; [self update:U_TOTAL]; return self; } /* read a MolViewer file */ -(BOOL)readMVfrag:(char *)name { FILE *in; int i; char s[80],*bpath; struct MVatm atm; struct MVbnd bnd; unsigned short n; bpath=(char *)[[NXBundle mainBundle] directory]; sprintf(s,"%s/%s.mvb",bpath,name); in=fopen(s,"r"); if (in==NULL) return NO; fread(&i,4,1,in); if (strncmp((char *)&i,"MVB",4)) { fclose(in); return NO; } if (frag.atom!=NULL) { free(frag.atom); free(frag.coord); free(frag.lb); } fread(&n,2,1,in); frag.numa=n; fread(&n,2,1,in); frag.numlb=n; fread(&i,4,1,in); /* for future expansion */ frag.atom=malloc(sizeof(Atom)*(frag.numa+1)); frag.coord=malloc(sizeof(RtPoint)*(frag.numa+1)); frag.lb=malloc(sizeof(struct LBOND)*frag.numlb); for (i=0; i<frag.numa; i++) { fread(&atm,18,1,in); if (atm.anum==255) atm.anum=84; frag.atom[i].anum=atm.anum; if (atm.type<32) atm.type=32; frag.atom[i].type=atm.type; frag.coord[i][0]=atm.x; frag.coord[i][1]=atm.y; frag.coord[i][2]=atm.z; frag.atom[i].charge=atm.c; frag.atom[i].numb=0; frag.atom[i].tag=frag.atom[i].sel=frag.atom[i].lab=0; } for (i=0; i<frag.numlb; i++) { fread(&bnd,5,1,in); frag.lb[i].n1=bnd.n1; frag.lb[i].n2=bnd.n2; frag.lb[i].type=bnd.type; frag.lb[i].d=frag.lb[i].w=0; } mol[nmol].ell=NULL; fclose(in); [self asnFragBnd]; return YES; } /* read a molviewer text file */ -readMVTX:(FILE *)infile { FILE *in; int i,j,k; short n1,n2; char s[100]; [self remMol:self]; in=infile; fscanf(in," "); s[0]='#'; while (s[0]=='#') fgets(s,100,in); sscanf(s," %hd %hd",&n1,&n2); mol[nmol].numa=n1; mol[nmol].numlb=n2; mol[nmol].atom=malloc(sizeof(Atom)*(mol[nmol].numa+1)); mol[nmol].coord=malloc(sizeof(RtPoint)*(mol[nmol].numa+1)); mol[nmol].lb=malloc(sizeof(struct LBOND)*mol[nmol].numlb); for (i=0; i<mol[nmol].numa; i++) { fgets(s,100,in); j=[self parseAtom:s :&mol[nmol].atom[i].anum :&mol[nmol].atom[i].type]; k=sscanf(&s[j]," %f %f %f %f %d",&mol[nmol].coord[i][0], &mol[nmol].coord[i][1],&mol[nmol].coord[i][2], &mol[nmol].atom[i].charge, &j); if (k!=5) mol[nmol].atom[i].lab=0; else mol[nmol].atom[i].lab=j; mol[nmol].atom[i].numb=0; mol[nmol].atom[i].tag=mol[nmol].atom[i].sel=0; } for (i=0; i<mol[nmol].numlb; i++) { fgets(s,100,in); sscanf(s," %d %d %c",&mol[nmol].lb[i].n1, &mol[nmol].lb[i].n2,&mol[nmol].lb[i].type); mol[nmol].lb[i].d=mol[nmol].lb[i].w=0.0; } if (fgets(s,100,in)!=NULL&&strncmp(s,"ELLIPSE",7)==0) { mol[nmol].ell=malloc(sizeof(Ellipse)*(mol[nmol].numa+1)); for (i=0; i<mol[nmol].numa; i++) { fgets(s,100,in); sscanf(s," %f %f %f %f %f %f",&mol[nmol].ell[i].x, &mol[nmol].ell[i].y,&mol[nmol].ell[i].z,&mol[nmol].ell[i].dz, &mol[nmol].ell[i].dx,&mol[nmol].ell[i].dz2); } } else mol[nmol].ell=NULL; fclose(in); [self asnBnd:nmol]; nmol++; [self centerMol:0]; [inspector newMol]; [self update:U_TOTAL]; return self; } /* read a Protein Data Bank molecule */ -readPdb:sender { id panel; FILE *in; int i,j,k,l,nb,b[11]; float x,y,z; char s[101],t[40]; RtPoint *a1,*a2; if (sender==nil) in=fopen(f2open,"r"); else { panel=[[OpenPanel new] allowMultipleFiles:NO]; if (![panel runModalForTypes:NULL]) return self; in=fopen([panel filename],"r"); } if (in==NULL) return self; /* get rid of the current molecule */ [self remMol:self]; i=nb=0; k=1; x=y=z=0; /* first pass, count atoms and bonds */ while (fgets(s,100,in)!=NULL) { /* still uses CONECT records to generate bonds. Needs to be replaced */ if (strncmp(s,"CONECT",6)==0) { s[70]=0; j=sscanf(&s[7]," %d %d %d %d %d %d %d %d %d %d",&k,&k,&k,&k,&k,&k,&k,&k,&k,&k); nb+=j-1; } if (k && (strncmp(s,"ATOM",4)==0 || strncmp(s,"HETATM",6)==0)) { sscanf(&s[7]," %d",&i); } /* only reads the first structure in the file, TER terminates reading */ if (strncmp(s,"TER",3)==0) k=0; } rewind(in); mol[nmol].numa=i; mol[nmol].atom=malloc(sizeof(Atom)*(i+1)); mol[nmol].coord=malloc(sizeof(RtPoint)*(i+1)); mol[nmol].numlb=nb; if (nb!=0) mol[nmol].lb=malloc(sizeof(struct LBOND)*(nb+1)); i=nb=0; k=0; /* ok, this time actually read the data */ while (fgets(s,100,in)!=NULL) { l=strlen(s); if (strncmp(s,"CONECT",6)==0) { s[70]=0; for (i=7,j=0; i<58; i+=5,j++) { s[i+4]=0; if (i>=l) { b[j]=-1; continue; } if (sscanf(&s[i]," %d",b+j)!=1) b[j]= -1; } a1=(RtPoint *)&mol[nmol].coord[b[0]-1][0]; for (i=1; i<11; i++) { if (nb>mol[nmol].numlb) { [self error:"Bad bond count" :NO]; continue; } if (b[i]==-1) continue; a2=(RtPoint *)&mol[nmol].coord[b[i]-1][0]; if (b[i]>mol[nmol].numa || b[i+1]>mol[nmol].numa) continue; for (j=0; j<nb; j++) { if (mol[nmol].lb[j].n1==b[0]-1&&mol[nmol].lb[j].n2==b[i]-1) break; if (mol[nmol].lb[j].n2==b[0]-1&&mol[nmol].lb[j].n1==b[i]-1) break; } if (j!=nb) continue; mol[nmol].lb[nb].n1=b[0]-1; mol[nmol].lb[nb].n2=b[i]-1; mol[nmol].lb[nb].w=1.0; mol[nmol].lb[nb].d=sqrt(SQR(a1[0]-a2[0])+SQR(a1[1]-a2[1])+ SQR(a1[2]-a2[2])); nb++; } continue; } if (k) continue; if (strncmp(s,"TER",3)==0) k=1; if (strncmp(s,"ATOM",4)!=0 && strncmp(s,"HETATM",6)!=0) continue; sscanf(&s[6]," %d",&i); i--; /* try to find this atom in the atomic table */ for (j=0; j<NEL; j++) if(s[12]==elinfo[j].name[0]&&s[13]==elinfo[j].name[1]) break; if (j==NEL) { for (j=0; j<NEL; j++) if(s[13]==elinfo[j].name[0] && elinfo[j].name[1]==' ') break; } if (j==NEL) { sprintf(t,"Can't find atom:%c%c",s[12],s[13]); [self error:t :NO]; } mol[0].atom[i].type=' '; if (s[14]=='A'&&j==5) mol[0].atom[i].type='T'; else if (s[14]==' '&&j==6) mol[0].atom[i].type='M'; mol[nmol].atom[i].anum=j; mol[nmol].atom[i].numb=0; mol[nmol].atom[i].tag=mol[nmol].atom[i].sel=mol[nmol].atom[i].lab=0; sscanf(&s[30]," %f %f %f",&mol[nmol].coord[i][0], &mol[nmol].coord[i][1],&mol[nmol].coord[i][2]); /* find the "average" position for centering */ x+=mol[nmol].coord[i][0]; y+=mol[nmol].coord[i][1]; z+=mol[nmol].coord[i][2]; } mol[nmol].numlb=nb; x/=(float)mol[nmol].numa; y/=(float)mol[nmol].numa; z/=(float)mol[nmol].numa; /* move origin to center of molecule */ for (i=0; i<mol[nmol].numa; i++) { mol[nmol].coord[i][0]-=x; mol[nmol].coord[i][1]-=y; mol[nmol].coord[i][2]-=z; } mol[nmol].ell=NULL; fclose(in); [self asnBnd:nmol]; nmol++; [inspector newMol]; [self update:U_TOTAL]; return self; } /* read an Alchemy format file */ -readAlch:sender { id panel; FILE *in; int i,j,k,l; float x,y,z,x1,y1,z1; float *a1,*a2; char c[4],s[101],t[40]; x=y=z=0.0; if (sender==nil) in=fopen(f2open,"r"); else { panel=[[OpenPanel new] allowMultipleFiles:NO]; if (![panel runModalForTypes:NULL]) return self; in=fopen([panel filename],"r"); } if (in==NULL) return self; [self remMol:self]; /* Ah Ha! Alchemy has a header with an atom/bond count */ fgets(s,100,in); sscanf(s," %d",&mol[nmol].numa); sscanf(&s[13]," %d",&mol[nmol].numlb); mol[nmol].atom=malloc(sizeof(Atom)*(mol[nmol].numa+1)); mol[nmol].coord=malloc(sizeof(RtPoint)*(mol[nmol].numa+1)); for (i=0; i<mol[nmol].numa; i++) { fgets(s,100,in); if (s[0]=='#') { i=0; continue; } sscanf(s," %d %c%c%c%c %f %f %f",&j,c,&c[1],&c[2],&c[3],&x1,&y1,&z1); j--; mol[nmol].coord[j][0]=x1; mol[nmol].coord[j][1]=y1; mol[nmol].coord[j][2]=z1; mol[nmol].atom[j].charge=0.0; mol[nmol].atom[j].numb=0; /* mol[nmol].atom[j].cnum=i; /* atom # within residue (for proteins) */ /* try to identify the atoms. */ if (c[0]=='~') c[0]='*'; for (k=0; k<NEL; k++) { if ((c[1]>='0'&&c[1]<='9')||c[1]==' '||c[2]!=' ') { if (c[0]==elinfo[k].name[0]&&elinfo[k].name[1]==' ') break; } else if (c[0]==elinfo[k].name[0]&&c[1]==elinfo[k].name[1]) break; } if (k==NEL) { sprintf(t,"Can't find atom:%c\n",c[0]); [self error:t :NO]; sleep(1); } mol[nmol].atom[j].anum=k; /* relabel atom types */ if (islower(c[1])) c[1]=toupper(c[1]); if (islower(c[2])) c[2]=toupper(c[2]); if (islower(c[3])) c[3]=toupper(c[3]); mol[nmol].atom[j].type=' '; if (c[1]=='A' && c[2]=='R') { if (k==5||k==6) mol[nmol].atom[j].type='A'; } else if (c[1]=='A' && c[2]=='M' && k==6) mol[nmol].atom[j].type='M'; else if (k==5 && c[1]=='3') mol[nmol].atom[j].type='T'; mol[nmol].atom[j].tag=mol[nmol].atom[j].sel=mol[nmol].atom[j].lab=0; x+=x1; /* keep running total for finding average x,y,z */ y+=y1; z+=z1; } j++; x/=(float)j; y/=(float)j; z/=(float)j; /* move origin to center of molecule */ for (i=0; i<mol[nmol].numa; i++) { mol[nmol].coord[i][0]-=x; mol[nmol].coord[i][1]-=y; mol[nmol].coord[i][2]-=z; } mol[nmol].lb=malloc(sizeof(struct LBOND)*mol[nmol].numlb); /* now read the bonds */ j=0; for (i=0; i<mol[nmol].numlb; i++) { fgets(s,100,in); sscanf(s,"%d %d %d",&j,&k,&l); mol[nmol].lb[i].n1=k-1; mol[nmol].lb[i].n2=l-1; a1=(float *)&mol[nmol].coord[k-1][0]; a2=(float *)&mol[nmol].coord[l-1][0]; mol[nmol].lb[i].w=1.0; mol[nmol].lb[i].type=s[19]; x=sqrt(SQR(a1[0]-a2[0])+SQR(a1[1]-a2[1])+SQR(a1[2]-a2[2])); if (x>2.2) x=1.8; mol[nmol].lb[i].d=x; } /*for (i=0; i<mol[0].numa; i++) { mol[0].coord[i][1]*=.5; mol[0].coord[i][0]=frand(-4.0,4.0); mol[0].coord[i][2]=frand(-4.0,4.0); }*/ mol[nmol].ell=NULL; [self asnBnd:nmol]; nmol++; fclose(in); [inspector newMol]; [self update:U_TOTAL]; return self; } /* delete the current molecule from memory */ -remMol:sender { if (nmol==0) return self; nmol--; if (mol[nmol].numa!=0) free(mol[nmol].atom); if (mol[nmol].numlb!=0) free(mol[nmol].lb); if (mol[nmol].ell!=NULL) free(mol[nmol].ell); if (mol[nmol].verts!=NULL) { free(mol[nmol].verts); mol[nmol].verts=NULL; } if (mol[nmol].nverts!=NULL) { free(mol[nmol].nverts); mol[nmol].nverts=NULL; } mol[nmol].numa=mol[nmol].numlb=0; [self update:U_TOTAL]; return self; } /* something has changed, so make sure all the displays are correct */ -update:(int)level { [nrgobj update:nmol level:level]; if (level>=U_BONDS) [self setupMol]; [molView setData:mol :nmol :elinfo]; [inspector update:nmol :level]; return self; } /* generates a data table to speed up drawing */ -setupMol { int i,j,k; for (i=0; i<nmol; i++) { k=mol[i].numlb; if (mol[i].verts!=NULL) free(mol[i].verts); if (mol[i].nverts!=NULL) free(mol[i].nverts); mol[i].verts=malloc(sizeof(RtInt)*2*(k+1)); mol[i].nverts=malloc(sizeof(RtInt)*(k+1)); for (j=0; j<k; j++) { mol[i].nverts[j]=2; mol[i].verts[j*2]=mol[i].lb[j].n1; mol[i].verts[j*2+1]=mol[i].lb[j].n2; } } return self; } /* Write an Alchemy format file */ -saveAlch:sender { int i,j; static id savePanel=nil; FILE *out; char s[120],type[3]; if (savePanel==nil) savePanel=[SavePanel new]; [savePanel setRequiredFileType:"mol"]; if (![savePanel runModal]) return self; out=fopen([savePanel filename],"w"); fprintf(out,"%5d ATOMS, %5d BONDS, 0 CHARGES, UNK\n", mol[0].numa,mol[0].numlb); for (i=0; i<mol[0].numa; i++) { if (elinfo[mol[0].atom[i].anum].name[1]==' ') sprintf(s,"%c",elinfo[mol[0].atom[i].anum].name[0]); else sprintf(s,"%c%c",elinfo[mol[0].atom[i].anum].name[0], elinfo[mol[0].atom[i].anum].name[1]); switch(mol[0].atom[i].anum) { case 0: type[0]=0; break; case 5: if (mol[0].atom[i].type!=' '&&mol[0].atom[i].type!='T') { type[0]='A'; type[1]='R'; type[2]=0; } else { type[0]=mol[0].atom[i].numb+'0'-1; type[1]=0; } break; case 6: if (mol[0].atom[i].type=='M') { type[0]='A'; type[1]='M'; type[2]=0; } else if (mol[0].atom[i].type=='A'|| mol[0].atom[i].type=='B'|| mol[0].atom[i].type=='C') { type[0]='A'; type[1]='R'; type[2]=0; } else { type[0]=mol[0].atom[i].numb+'0'-1; type[1]=0; } break; case 7: case 15: if (mol[0].atom[i].numb==1) { type[0]='2'; type[1]=0; } else { type[0]='3'; type[1]=0; } break; case 14: if (mol[0].atom[i].numb==5) { type[0]='4'; type[1]=0; } else { type[0]='3'; type[1]=0; } break; case 85: type[0]='P'; type[1]=0; break; default: type[0]=0; } strcat(s,type); for (j=strlen(type); j<4; j++) strcat(s," "); fprintf(out,"%5d %s%8.4f %8.4f %8.4f 0.0000\n", i+1,s,mol[0].coord[i][0],mol[0].coord[i][1], mol[0].coord[i][2]); } for (i=0; i<mol[0].numlb; i++) { if (mol[0].lb[i].type=='S') strcpy(s,"SINGLE"); else strcpy(s,"DOUBLE"); fprintf(out,"%5d %5d %5d %s\n",i+1,mol[0].lb[i].n1+1,mol[0].lb[i].n2+1,s); } fclose(out); return self; } /* Write an MV binary format file */ -saveMVB:sender { int i; static id savePanel=nil; FILE *out; struct MVatm atm; struct MVbnd bnd; unsigned short n; if (savePanel==nil) savePanel=[SavePanel new]; [savePanel setRequiredFileType:"mvb"]; if (![savePanel runModal]) return self; out=fopen([savePanel filename],"w"); /* magic # */ fwrite("MVB",4,1,out); /* # atoms and bonds */ n=mol[0].numa; fwrite(&n,sizeof(short),1,out); n=mol[0].numlb; fwrite(&n,sizeof(short),1,out); /* the next 2 are for future expansion */ n=0; fwrite(&n,sizeof(short),1,out); fwrite(&n,sizeof(short),1,out); for (i=0; i<mol[0].numa; i++) { atm.anum=mol[0].atom[i].anum; if (atm.anum==84) atm.anum=255; atm.type=mol[0].atom[i].type; atm.x=mol[0].coord[i][0]; atm.y=mol[0].coord[i][1]; atm.z=mol[0].coord[i][2]; atm.c=mol[0].atom[i].charge; fwrite(&atm,18,1,out); } for (i=0; i<mol[0].numlb; i++) { bnd.n1=mol[0].lb[i].n1; bnd.n2=mol[0].lb[i].n2; bnd.type=mol[0].lb[i].type; fwrite(&bnd,5,1,out); } fclose(out); return self; } /* Write an MV text format file */ -saveMVTX:sender { int i; static id savePanel=nil; FILE *out; unsigned short n; if (savePanel==nil) savePanel=[SavePanel new]; [savePanel setRequiredFileType:"mvt"]; if (![savePanel runModal]) return self; out=fopen([savePanel filename],"w"); /* # atoms and bonds */ fprintf(out,"MVTX\n%d %d\n",mol[0].numa,mol[0].numlb); for (i=0; i<mol[0].numa; i++) { n=mol[0].atom[i].anum; if (elinfo[n].name[1]==' ') fprintf(out,"%c%c\t",elinfo[n].name[0], mol[0].atom[i].type); else fprintf(out,"%c%c%c\t",elinfo[n].name[0],elinfo[n].name[1], mol[0].atom[i].type); fprintf(out,"%f\t%f\t%f\t%f\t%1d\n",mol[0].coord[i][0], mol[0].coord[i][1],mol[0].coord[i][2],mol[0].atom[i].charge, mol[0].atom[i].lab); } for (i=0; i<mol[0].numlb; i++) { fprintf(out,"%d\t%d\t%c\n",mol[0].lb[i].n1,mol[0].lb[i].n2,mol[0].lb[i].type); } if (mol[0].ell!=NULL) { fprintf(out,"ELLIPSE\n"); for (i=0; i<mol[0].numa; i++) { fprintf(out,"%f\t%f\t%f\t%f\t%f\t%f\n",mol[0].ell[i].x, mol[0].ell[i].y,mol[0].ell[i].z,mol[0].ell[i].dz, mol[0].ell[i].dx,mol[0].ell[i].dz2); } } fclose(out); return self; } /* attach a loaded fragment to the current molecule. n=* atom in mol[0] */ /* nf=* atom in frag. len is the bond length to use */ -attachFrag:(int)n :(int)nf :(float)ln :(char)type { int nc,ncf,s1,s2,i; float v1[3],v2[3],v3[3]; float ang,t; nc=mol[0].atom[n].bnd[0]; ncf=frag.atom[nf].bnd[0]; /* find the 2 bond vectors */ sub(&mol[0].coord[n][0],&mol[0].coord[nc][0],v1); sub(&frag.coord[ncf][0],&frag.coord[nf][0],v2); /* find a vector to rotate frag so the bonds point in the same dir */ cross(v1,v2,v3); /*printf("(%f,%f,%f) (%f,%f,%f) -> ",v1[0],v1[1],v1[2],v2[0],v2[1],v2[2]);*/ ang=asin(len(v3)/(len(v1)*len(v2))); if (dot(v1,v2)<0.0) { if (ang>0.0) ang=M_PI-ang; else ang= -M_PI-ang; } [self rotateFragAboutV:ang :v3]; sub(&frag.coord[ncf][0],&frag.coord[nf][0],v2); /*printf("(%f,%f,%f)\n",v2[0],v2[1],v2[2]); printf("%f %f\n",ang/DRC,acos(dot(v1,v2)/(len(v1)*len(v2)))/DRC);*/ t=ln/len(v1); mul(v1,t); add2(v1,mol[0].coord[nc]); sub2(v1,frag.coord[ncf]); s1=mol[0].numa; s2=mol[0].numlb; mol[0].numa+=frag.numa; mol[0].numlb+=frag.numlb; mol[0].atom=realloc(mol[0].atom,sizeof(Atom)*(mol[0].numa+1)); mol[0].coord=realloc(mol[0].coord,sizeof(RtPoint)*(mol[0].numa+1)); mol[0].lb=realloc(mol[0].lb,sizeof(struct LBOND)*(mol[0].numlb+1)); for (i=0; i<frag.numa; i++) { mol[0].atom[i+s1]=frag.atom[i]; add(v1,frag.coord[i],mol[0].coord[i+s1]); } for (i=0; i<frag.numlb; i++) { mol[0].lb[i+s2].n1=frag.lb[i].n1+s1; mol[0].lb[i+s2].n2=frag.lb[i].n2+s1; mol[0].lb[i+s2].type=frag.lb[i].type; } /* handled in delatom below [self delBond:n :nc]; [self delBond:nf+s1 :ncf+s1]; */ mol[0].lb[mol[0].numlb].n1=nc; mol[0].lb[mol[0].numlb].n2=ncf+s1; mol[0].lb[mol[0].numlb].type=type; mol[0].numlb++; [self delAtom:n]; [self delAtom:nf+s1-1]; [self asnBnd:0]; return self; } /* deletes an atom from mol[0]. You MUST do an asnBnd sometime after this. */ -delAtom:(int)n { int i,j; mol[0].numa--; for (i=n; i<mol[0].numa; i++) { mol[0].atom[i]=mol[0].atom[i+1]; mol[0].coord[i][0]=mol[0].coord[i+1][0]; mol[0].coord[i][1]=mol[0].coord[i+1][1]; mol[0].coord[i][2]=mol[0].coord[i+1][2]; } j=0; for (i=0; i<mol[0].numlb; i++) { if (mol[0].lb[i].n1==n || mol[0].lb[i].n2==n) continue; if (mol[0].lb[i].n1>n) mol[0].lb[i].n1--; if (mol[0].lb[i].n2>n) mol[0].lb[i].n2--; if (i!=j) mol[0].lb[j]=mol[0].lb[i]; j++; } mol[0].numlb=j; return self; } -makeProt:(int)n :(int *)aa { int i,j,k,fn=0,fc=0,nm=0,nm2=0,fn2=0,fc2=0; [self remMol:self]; /* read 1st acid & copy -> mol[0] */ [self readMVfrag:acid[aa[0]].d2]; mol[0].atom=malloc(sizeof(Atom)*(frag.numa+1)); mol[0].coord=malloc(sizeof(RtPoint)*(frag.numa+1)); mol[0].lb=malloc(sizeof(struct LBOND)*frag.numlb); mol[0].numlb=frag.numlb; mol[0].numa=frag.numa; for (i=0; i<frag.numa; i++) { mol[0].atom[i]=frag.atom[i]; mol[0].coord[i][0]=frag.coord[i][0]; mol[0].coord[i][1]=frag.coord[i][1]; mol[0].coord[i][2]=frag.coord[i][2]; } for (i=0; i<frag.numlb; i++) mol[0].lb[i]=frag.lb[i]; for (j=1; j<n; j++) { if ([self readMVfrag:acid[aa[j]].d2]==NO) continue; /* locate the free bonds in mol[0]*/ for (i=0; i<mol[0].numa; i++) { if (mol[0].atom[i].anum==84) { k=mol[0].atom[i].bnd[0]; if (mol[0].atom[k].anum==6&&mol[0].atom[k].type=='M') fn=i; else fc=i; } } /* locate the free bonds in frag */ for (nm2=0; nm2<frag.numa; nm2++) if (frag.atom[nm2].anum==6&&frag.atom[nm2].type=='M') break; for (i=0; i<frag.numa; i++) { if (frag.atom[i].anum==84) { if (frag.atom[i].bnd[0]==nm2) fn2=i; else fc2=i; } } [self attachFrag:fc :fn2 :1.335 :'M']; } for (i=0; i<mol[0].numa; i++) { if (mol[0].atom[i].anum==84) [self delAtom:i]; } mol[nmol].ell=NULL; [self asnBnd:0]; nmol=1; [self centerMol:0]; [self update:U_TOTAL]; return self; } /* copy bond table to redundant bond information in Atom structure */ -asnBnd:(int)nm { int i,j; int nlb,n1,n2; float l; for(i=0; i<mol[nm].numa; i++) mol[nm].atom[i].numb=0; nlb=mol[nm].numlb; for (i=0; i<nlb; i++) { n1=mol[nm].lb[i].n1; n2=mol[nm].lb[i].n2; l=mol[nm].lb[i].d; j=mol[nm].atom[n1].numb++; if (j>3) { j=3; /*printf("Too many bonds\n");*/ } mol[nm].atom[n1].bnd[j]=n2; /* mol[nm].atom[n1].bndl[j]=l; */ j=mol[nm].atom[n2].numb++; if (j>3) { j=3; /*printf("Too many bonds\n");*/ } mol[nm].atom[n2].bnd[j]=n1; /* mol[nm].atom[n2].bndl[j]=l; */ } for (i=0; i<mol[nm].numa; i++) { for (j=mol[nm].atom[i].numb; j<4; j++) mol[nm].atom[i].bnd[j]=-1; } if (nm==0) chi=[self calcChi:self]; return self; } /* copy bond table to redundant bond information in Atom structure */ -asnFragBnd { int i,j; int lb,n1,n2; float l; lb=frag.numlb; for (i=0; i<lb; i++) { n1=frag.lb[i].n1; n2=frag.lb[i].n2; l=frag.lb[i].d; j=frag.atom[n1].numb++; if (j>3) { j=3; /*printf("Too many bonds\n");*/ } frag.atom[n1].bnd[j]=n2; /* frag.atom[n1].bndl[j]=l; */ j=frag.atom[n2].numb++; if (j>3) { j=3; /*printf("Too many bonds\n");*/ } frag.atom[n2].bnd[j]=n1; /* frag.atom[n2].bndl[j]=l; */ } for (i=0; i<frag.numa; i++) { for (j=frag.atom[i].numb; j<4; j++) frag.atom[i].bnd[j]=-1; } return self; } /* locate and identify amino acids in the current molecule */ -findAA:sender { int h,i,j,k,l,m,n,na; int cnt[16]; char s[80]; struct AMINO tmp; /* This is difficult to follow, but you'll just have to wade through it */ /* if you need to change something. */ na=0; for (i=0; i<mol[0].numa; i++) { if (mol[0].atom[i].anum!=6 || mol[0].atom[i].type!='M') continue; mol[0].am[na].n=i; /* assign N */ mol[0].atom[i].tag=1; mol[0].am[na].cn= -1; mol[0].am[na].h= -1; /* assign h,c1, and cn */ m= -1; for (j=0; j<3; j++) { k=mol[0].atom[i].bnd[j]; if (k==-1) continue; /* h is attached to n */ if (mol[0].atom[k].anum==0) mol[0].am[na].h=k; /* cn has an O attached, c1 has an H and the sidechain */ if (mol[0].atom[k].anum==5) { if (mol[0].atom[k].type=='T') { for (l=h=0; l<mol[0].atom[k].numb; l++) { n=mol[0].atom[k].bnd[l]; if (mol[0].atom[n].anum==0) h++; } if (h==1) mol[0].am[na].c1=k; else m=k; } else if (mol[0].atom[k].type==' ') mol[0].am[na].cn=k; else [self error:"Bad carbon on Nam" :NO]; } } if (m!=-1) { if (mol[0].am[na].h==-1) mol[0].am[na].h=m; else mol[0].am[na].c1=m; } /* assign cn2 */ k=mol[0].am[na].cn; if (k==-1) mol[0].am[na].cn2= -1; else { /* cn2 is attached to N,H,and a sc */ for (j=0; j<mol[0].atom[k].numb; j++) { n=mol[0].atom[k].bnd[j]; if (mol[0].atom[n].anum==5 && mol[0].atom[n].type=='T') mol[0].am[na].cn2=n; } } /* assign c2 */ k=mol[0].am[na].c1; mol[0].atom[k].tag=1; /* c2 is attached to O */ for (j=0; j<mol[0].atom[k].numb; j++) { n=mol[0].atom[k].bnd[j]; if (mol[0].atom[n].anum==5 && mol[0].atom[n].type==' ') mol[0].am[na].c2=n; } /* assign o */ k=mol[0].am[na].c2; for (j=0; j<mol[0].atom[k].numb; j++) { n=mol[0].atom[k].bnd[j]; if (mol[0].atom[n].anum==7) mol[0].am[na].o=n; } /* assign sc */ k=mol[0].am[na].c1; h=0; for (j=0; j<4; j++) { n=mol[0].atom[k].bnd[j]; if (n==mol[0].am[na].c2||n==mol[0].am[na].cn) continue; if (mol[0].atom[n].anum==0) { h=n; continue; } if (mol[0].atom[n].anum==5) { mol[0].am[na].sc=n; break; } } /* if glycine, an h is chosen randomly */ if (j==4) mol[0].am[na].sc=h; /* assign sc2 */ k=mol[0].am[na].sc; for (j=0; j<4; j++) { n=mol[0].atom[k].bnd[j]; if (n==mol[0].am[na].c1||n==-1) continue; if (mol[0].atom[n].anum==0) { h=n; continue; } else { mol[0].am[na].sc2=n; break; } } if (j==4) mol[0].am[na].sc2=h; na++; if (na==MAXAA) { na--; [self error:"Max # of amino acids exceeded" :NO]; } } mol[0].numaa=na; for (i=0; i<mol[0].numaa; i++) { /* calculate the dihedrals */ mol[0].am[i].sel=0; 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]; 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; 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; mol[0].am[i].chi=[self dihedral:mol[0].am[i].n :mol[0].am[i].c1 :mol[0].am[i].sc :mol[0].am[i].sc2]; mol[0].am[i].anum=0; if (mol[0].am[i].psi< -M_PI) mol[0].am[i].psi+=M_PI; if (mol[0].am[i].phi< -M_PI) mol[0].am[i].phi+=M_PI; /* try to identify the individual amino acids by counting and a little */ /* chain testing when the result is ambigous. Exotic amino acids */ /* might be incorrectly assigned. Proline is also incorrectly assigned. */ for (j=0; j<16; j++) cnt[j]=0; j=mol[0].am[i].sc; k=mol[0].am[i].h; if (mol[0].atom[j].anum==0) mol[0].am[i].anum=9; /* Gly */ else if (mol[0].atom[k].anum==5) mol[0].am[i].anum=5; /* Pro */ else { count(mol,cnt,mol[0].am[i].sc,0,15); for (k=0; k<AADEF; k++) { if (acid[k].nc==cnt[5]&&acid[k].ns==cnt[15]&&acid[k].no==cnt[7]&& acid[k].nn==cnt[6]) break; } if (k==AADEF) { mol[0].am[i].anum=0; sprintf(s,"Unknown amino acid: C=%d S=%d O=%d N=%d",cnt[5],cnt[15],cnt[7],cnt[6]); [self error:s :NO]; } else mol[0].am[i].anum=k; } if (mol[0].am[i].anum==3) { l=mol[0].am[i].sc; cnt[5]=0; for (k=0; k<mol[0].atom[l].numb; k++) if (mol[0].atom[mol[0].atom[l].bnd[k]].anum==5) cnt[5]++; if (cnt[5]==3) mol[0].am[i].anum=4; } } for (i=0; i<mol[0].numa; i++) mol[0].atom[i].tag=0; /* sort amino acids so N terminus is FIRST */ for (i=1; i<mol[0].numaa; i++) if (mol[0].am[i].cn==-1) swapa(0,i); for (i=1; i<mol[0].numaa; i++) { for (j=i+2; j<mol[0].numaa; j++) { if (mol[0].am[i-1].c2==mol[0].am[j].cn) swapa(i,j); } } return self; } /* recursive routine to decend the sidechain and count the number of */ /* atoms of each type */ void count (Molecule *mol,int *cnt, int n1,short lev,short maxlev) { int i,j; if (lev>maxlev||mol[0].atom[n1].tag) return; j=mol[0].atom[n1].anum; if (j<16) cnt[j]++; mol[0].atom[n1].tag=1; for (i=0; i<mol[0].atom[n1].numb; i++) { count(mol,cnt,mol[0].atom[n1].bnd[i],lev+1,maxlev); } return; } /* test routine (not used) */ -test { int i; for (i=0; i<mol[0].numa; i++) printf("%c%c%c,",elinfo[mol[0].atom[i].anum].name[0],elinfo[mol[0].atom[i].anum].name[1],mol[0].atom[i].type); printf(" %d\n",mol[0].numa); for (i=0; i<mol[0].numlb; i++) printf("%d-%d ",mol[0].lb[i].n1,mol[0].lb[i].n2); printf("%d\n",mol[0].numlb); for (i=0; i<frag.numa; i++) printf("%c%c%c,",elinfo[frag.atom[i].anum].name[0],elinfo[frag.atom[i].anum].name[1],frag.atom[i].type); printf(" %d\n",frag.numa); for (i=0; i<frag.numlb; i++) printf("%d-%d ",frag.lb[i].n1,frag.lb[i].n2); printf("%d\n",frag.numlb); return self; } /* set amino acid dihedral angles */ -setOmega:(float)ang { int i; float a,b; a=ang*DRC; for (i=0; i<mol[0].numaa; i++) { if (mol[0].am[i].cn==-1||(!mol[0].am[i].sel)) continue; mol[0].am[i].omega=a; b=[self dihedral:mol[0].am[i].c1 :mol[0].am[i].n :mol[0].am[i].cn :mol[0].am[i].cn2]; [self rotateAbout:a-b :mol[0].am[i].c1 :mol[0].am[i].n :mol[0].am[i].cn]; } if (updFlag) [self centerMol:0]; return self; } -setPsi:(float)ang { int i; float a,b; a=ang+180.0; if (a>180.0) a-=360.0; a*=DRC; for (i=0; i<mol[0].numaa; i++) { if (!mol[0].am[i].sel) continue; mol[0].am[i].psi=ang*DRC; b=[self dihedral:mol[0].am[i].n :mol[0].am[i].c1 :mol[0].am[i].c2 :mol[0].am[i].o]; [self rotateAbout:a-b :mol[0].am[i].n :mol[0].am[i].c1 :mol[0].am[i].c2]; } if (updFlag) [self centerMol:0]; return self; } -setPhi:(float)ang { int i; float a,b; a=ang+180.0; if (a>180.0) a-=360.0; a*=DRC; for (i=0; i<mol[0].numaa; i++) { if (!mol[0].am[i].sel) continue; mol[0].am[i].phi=ang*DRC; b=[self dihedral:mol[0].am[i].c2 :mol[0].am[i].c1 :mol[0].am[i].n :mol[0].am[i].h]; [self rotateAbout:a-b :mol[0].am[i].c2 :mol[0].am[i].c1 :mol[0].am[i].n]; } if (updFlag) [self centerMol:0]; return self; } -setChi:(float)ang { int i; float a,b; a=ang; if (a>180.0) a-=360.0; a*=DRC; for (i=0; i<mol[0].numaa; i++) { if (!mol[0].am[i].sel) continue; if (mol[0].am[i].sc2==-1) continue; mol[0].am[i].chi=ang*DRC; b=[self dihedral:mol[0].am[i].n :mol[0].am[i].c1 :mol[0].am[i].sc :mol[0].am[i].sc2]; [self rotateAbout:a-b :mol[0].am[i].n :mol[0].am[i].c1 :mol[0].am[i].sc]; } if (updFlag) [self centerMol:0]; return self; } /* take 4 atoms and calculate a dihedral angle. The math is right, but */ /* it would take too much space to explain it. */ -(float)dihedral:(int)n1 :(int)n2 :(int)n3 :(int)n4 { float cb[3],ca[3],bd[3],bc[3],v1[3],v2[3],v3[3]; float t,t2; if (n1==-1||n2==-1||n3==-1||n4==-1) return 0.0; sub(&mol[0].coord[n2][0],&mol[0].coord[n3][0],cb); sub(&mol[0].coord[n3][0],&mol[0].coord[n2][0],bc); sub(&mol[0].coord[n1][0],&mol[0].coord[n3][0],ca); sub(&mol[0].coord[n4][0],&mol[0].coord[n2][0],bd); cross(cb,ca,v1); cross(bd,bc,v2); cross(v1,v2,v3); t=dot(v3,bc)/len(bc); t2=dot(v1,v2); return(atan2(t,t2)); } /* calculate a bond angle */ -(float)bondAng:(int)n1 :(int)n2 :(int)n3; { float v1[3],v2[3],r; sub(&mol[0].coord[n1][0],&mol[0].coord[n2][0],v1); sub(&mol[0].coord[n3][0],&mol[0].coord[n2][0],v2); r=dot(v1,v2)/(len(v1)*len(v2)); return (acos(r)); } /* test if 2 atoms are bonded */ -(BOOL)bonded:(short)n1 :(short)n2 { int i; for (i=0; i<mol[0].atom[n1].numb; i++) if (mol[0].atom[n1].bnd[i]==n2) break; if (i==mol[0].atom[n1].numb) return(NO); return(YES); } /* center the molecule on the screen */ - centerMol:(int)m { float x,y,z,mx[9],a; int i,j,k,l; x=y=z=0.0; for (i=0; i<mol[m].numa; i++) { x+=mol[m].coord[i][0]; y+=mol[m].coord[i][1]; z+=mol[m].coord[i][2]; } x/=(float)mol[m].numa; y/=(float)mol[m].numa; z/=(float)mol[m].numa; for (i=0; i<mol[m].numa; i++) { mol[m].coord[i][0]-=x; mol[m].coord[i][1]-=y; mol[m].coord[i][2]-=z; } /*i=mol[m].am[4].o; j=mol[m].am[7].h; x=sqrt(SQR(mol[m].coord[i][0]-mol[m].coord[j][0])+ SQR(mol[m].coord[i][1]-mol[m].coord[j][1])+ SQR(mol[m].coord[i][2]-mol[m].coord[j][2])); i=mol[m].am[4].o; j=mol[m].am[6].h; y=sqrt(SQR(mol[m].coord[i][0]-mol[m].coord[j][0])+ SQR(mol[m].coord[i][1]-mol[m].coord[j][1])+ SQR(mol[m].coord[i][2]-mol[m].coord[j][2])); [hba setFloatValue:x]; [hb3 setFloatValue:y]; */ /* if there are >4 amino acids, align along Z */ if (mol[m].numaa>4) { if (mol[m].ell!=NULL) { free(mol[0].ell); mol[0].ell=NULL; [self error:"Ellipticity info lost!"]; } j=mol[m].am[1].n; k=mol[m].am[2].n; l=mol[m].am[3].n; y=(mol[m].coord[j][1]+mol[m].coord[k][1]+mol[m].coord[l][1])/3.0; z=(mol[m].coord[j][2]+mol[m].coord[k][2]+mol[m].coord[l][2])/3.0; a=atan2(z,y); mx[4]=mx[8]=cos(a); mx[5]=sin(a); mx[7]=-mx[5]; for (i=0; i<mol[m].numa; i++) { y=mol[m].coord[i][1]; z=mol[m].coord[i][2]; mol[m].coord[i][1]=y*mx[4]+z*mx[5]; mol[m].coord[i][2]=y*mx[7]+z*mx[8]; } y=(mol[m].coord[j][1]+mol[m].coord[k][1]+mol[m].coord[l][1])/3.0; x=(mol[m].coord[j][0]+mol[m].coord[k][0]+mol[m].coord[l][0])/3.0; a=atan2(x,y); mx[0]=mx[4]=cos(a); mx[1]=-sin(a); mx[3]=sin(a); for (i=0; i<mol[m].numa; i++) { y=mol[m].coord[i][1]; x=mol[m].coord[i][0]; mol[m].coord[i][0]=x*mx[0]+y*mx[1]; mol[m].coord[i][1]=x*mx[3]+y*mx[4]; } } [self update:U_POS]; return self; } /* aligns molecule along line connecting center and avg of sel. atoms */ -alignMol:sender { int i,j; float x,y,z,a,mx[9]; if (mol[0].ell!=NULL) { free(mol[0].ell); mol[0].ell=NULL; [self error:"Ellipticity info lost!" :NO]; } j=0; x=y=z=0.0; for (i=0; i<mol[0].numa; i++) { if (!mol[0].atom[i].sel) continue; y+=mol[0].coord[i][1]; z+=mol[0].coord[i][2]; j++; } if (j==0) { [self error:"Select atom(s) first." :NO]; return self; } y/=(float)j; z/=(float)j; a=atan2(z,y); mx[4]=mx[8]=cos(a); mx[5]=sin(a); mx[7]=-mx[5]; for (i=0; i<mol[0].numa; i++) { y=mol[0].coord[i][1]; z=mol[0].coord[i][2]; mol[0].coord[i][1]=y*mx[4]+z*mx[5]; mol[0].coord[i][2]=y*mx[7]+z*mx[8]; } j=0; x=y=z=0.0; for (i=0; i<mol[0].numa; i++) { if (!mol[0].atom[i].sel) continue; y+=mol[0].coord[i][1]; x+=mol[0].coord[i][0]; j++; } y/=(float)j; x/=(float)j; a=atan2(x,y); mx[0]=mx[4]=cos(a); mx[1]=-sin(a); mx[3]=sin(a); for (i=0; i<mol[0].numa; i++) { y=mol[0].coord[i][1]; x=mol[0].coord[i][0]; mol[0].coord[i][0]=x*mx[0]+y*mx[1]; mol[0].coord[i][1]=x*mx[3]+y*mx[4]; } [self update:U_POS]; return self; } /* set the view window to a variety of convenient sizes */ -setWinSize:sender { [vWindow sizeWindow:[[sender cellAt:0 :0] floatValue]+4.0 :[[sender cellAt:0 :1] floatValue]+2.0]; return self; } /* do a dihedral rotation about the line connecting n1 and n2 */ -rotateAbout:(float)a :(int)n0 :(int)n1 :(int)n2 { Atom *atom; int i; RtPoint n,r; float mx[9],m; a= -a; atom=mol[0].atom; for (i=0; i<mol[0].numa; i++) atom[i].tag=0; atom[n1].tag=1; tagatoms(atom,n2); atom[n1].tag=0; atom[n2].tag=0; atom[n0].tag=0; n[0]=mol[0].coord[n2][0]-mol[0].coord[n1][0]; /* n==vec from n1 to n2 */ n[1]=mol[0].coord[n2][1]-mol[0].coord[n1][1]; n[2]=mol[0].coord[n2][2]-mol[0].coord[n1][2]; m=sqrt(SQR(n[0])+SQR(n[1])+SQR(n[2])); n[0]/=m; /* n is unit vector */ n[1]/=m; n[2]/=m; m=1.0-cos(a); /* set up mx transformation matrix */ mx[0]=cos(a)+n[0]*n[0]*m; mx[1]=n[0]*n[1]*m+n[2]*sin(a); mx[2]=n[0]*n[2]*m-n[1]*sin(a); mx[3]=n[0]*n[1]*m-n[2]*sin(a); mx[4]=cos(a)+n[1]*n[1]*m; mx[5]=n[2]*n[1]*m+n[0]*sin(a); mx[6]=n[0]*n[2]*m+n[1]*sin(a); mx[7]=n[1]*n[2]*m-n[0]*sin(a); mx[8]=cos(a)+n[2]*n[2]*m; for (i=0; i<mol[0].numa; i++) { if (atom[i].tag==0) continue; r[0]=mol[0].coord[i][0]-mol[0].coord[n1][0]; r[1]=mol[0].coord[i][1]-mol[0].coord[n1][1]; r[2]=mol[0].coord[i][2]-mol[0].coord[n1][2]; mol[0].coord[i][0]=r[0]*mx[0]+r[1]*mx[1]+r[2]*mx[2]+mol[0].coord[n1][0]; mol[0].coord[i][1]=r[0]*mx[3]+r[1]*mx[4]+r[2]*mx[5]+mol[0].coord[n1][1]; mol[0].coord[i][2]=r[0]*mx[6]+r[1]*mx[7]+r[2]*mx[8]+mol[0].coord[n1][2]; atom[i].tag=0; } return self; } /* rotate the fragment about the vector n */ -rotateFragAboutV:(float)a :(float *)n { int i; float r[3]; float mx[9],m; m=sqrt(SQR(n[0])+SQR(n[1])+SQR(n[2])); n[0]/=m; /* n is unit vector */ n[1]/=m; n[2]/=m; m=1.0-cos(a); /* set up mx transformation matrix */ mx[0]=cos(a)+n[0]*n[0]*m; mx[1]=n[0]*n[1]*m+n[2]*sin(a); mx[2]=n[0]*n[2]*m-n[1]*sin(a); mx[3]=n[0]*n[1]*m-n[2]*sin(a); mx[4]=cos(a)+n[1]*n[1]*m; mx[5]=n[2]*n[1]*m+n[0]*sin(a); mx[6]=n[0]*n[2]*m+n[1]*sin(a); mx[7]=n[1]*n[2]*m-n[0]*sin(a); mx[8]=cos(a)+n[2]*n[2]*m; for (i=0; i<frag.numa; i++) { r[0]=frag.coord[i][0]; r[1]=frag.coord[i][1]; r[2]=frag.coord[i][2]; frag.coord[i][0]=r[0]*mx[0]+r[1]*mx[1]+r[2]*mx[2]; frag.coord[i][1]=r[0]*mx[3]+r[1]*mx[4]+r[2]*mx[5]; frag.coord[i][2]=r[0]*mx[6]+r[1]*mx[7]+r[2]*mx[8]; } return self; } -modBondLen:(float)l :(int)n1 :(int)n2 { Atom *atom; int i; float n[3]; float m; atom=mol[0].atom; for (i=0; i<mol[0].numa; i++) atom[i].tag=0; atom[n1].tag=1; tagatoms(atom,n2); atom[n1].tag=0; sub(&mol[0].coord[n1][0],&mol[0].coord[n2][0],n); m=len(n); n[0]*=(1.0-l/m); /* n is translation vector */ n[1]*=(1.0-l/m); n[2]*=(1.0-l/m); for (i=0; i<mol[0].numa; i++) { if (atom[i].tag==0) continue; mol[0].coord[i][0]+=n[0]; mol[0].coord[i][1]+=n[1]; mol[0].coord[i][2]+=n[2]; atom[i].tag=0; } return self; } /* rotate about the line (n1-n2)X(n3-n2) */ -modBondAng:(float)a :(int)n1 :(int)n2 :(int)n3 { Atom *atom; int i; float n[3],r[3],v1[3],v2[3]; float mx[9],m; a= -a; atom=mol[0].atom; for (i=0; i<mol[0].numa; i++) atom[i].tag=0; atom[n2].tag=1; tagatoms(atom,n3); atom[n2].tag=0; atom[n1].tag=0; sub(&mol[0].coord[n1][0],&mol[0].coord[n2][0],v1); sub(&mol[0].coord[n3][0],&mol[0].coord[n2][0],v2); cross(v1,v2,n); m=len(n); n[0]/=m; /* n is unit vector */ n[1]/=m; n[2]/=m; m=1.0-cos(a); /* set up mx transformation matrix */ mx[0]=cos(a)+n[0]*n[0]*m; mx[1]=n[0]*n[1]*m+n[2]*sin(a); mx[2]=n[0]*n[2]*m-n[1]*sin(a); mx[3]=n[0]*n[1]*m-n[2]*sin(a); mx[4]=cos(a)+n[1]*n[1]*m; mx[5]=n[2]*n[1]*m+n[0]*sin(a); mx[6]=n[0]*n[2]*m+n[1]*sin(a); mx[7]=n[1]*n[2]*m-n[0]*sin(a); mx[8]=cos(a)+n[2]*n[2]*m; for (i=0; i<mol[0].numa; i++) { if (atom[i].tag==0) continue; r[0]=mol[0].coord[i][0]-mol[0].coord[n2][0]; r[1]=mol[0].coord[i][1]-mol[0].coord[n2][1]; r[2]=mol[0].coord[i][2]-mol[0].coord[n2][2]; mol[0].coord[i][0]=r[0]*mx[0]+r[1]*mx[1]+r[2]*mx[2]+mol[0].coord[n2][0]; mol[0].coord[i][1]=r[0]*mx[3]+r[1]*mx[4]+r[2]*mx[5]+mol[0].coord[n2][1]; mol[0].coord[i][2]=r[0]*mx[6]+r[1]*mx[7]+r[2]*mx[8]+mol[0].coord[n2][2]; atom[i].tag=0; } return self; } /* select atoms from SelectView data Structure */ -setSelAtom:(struct SELDAT *)list { int i; for (i=0; i<mol[0].numa; i++) { if (list[i].sel) mol[0].atom[i].sel=1; else mol[0].atom[i].sel=0; } [molView display]; return self; } -setUpdFlag:sender { updFlag=[sender intValue]; return self; } /* takes a string with an atom name/type and gives anum/type */ -(int)parseAtom:(char *)s :(char *)anum :(char *)type { int j; char t[50]; if (isspace(s[1])) { type[0]=' '; if (s[0]=='~') s[0]='*'; for (j=0; j<NEL; j++) if(s[0]==elinfo[j].name[0]&&isspace(elinfo[j].name[1])) break; if (j==NEL) { sprintf(t,"Unidentified atom type %c.\n",s[0]); [self error:t :NO]; } anum[0]=j; return(2); } if (isspace(s[2])&&!islower(s[1])) { type[0]=s[1]; if (s[0]=='~') s[0]='*'; for (j=0; j<NEL; j++) if(s[0]==elinfo[j].name[0]&&isspace(elinfo[j].name[1])) break; if (j==NEL) { sprintf(t,"Unidentified atom type %c.\n",s[0]); [self error:t :NO]; } anum[0]=j; return(3); } if (isspace(s[2])) type[0]=' '; else type[0]=s[2]; for (j=0; j<NEL; j++) if(s[0]==elinfo[j].name[0]&&s[1]==elinfo[j].name[1]) break; if (j==NEL) { sprintf(t,"Unidentified atom type %c%c.\n",s[0],s[1]); [self error:t :NO]; } anum[0]=j; return(4); } /* recursive atom tagger for rotation routines */ void tagatoms(Atom *atom,int n) { int i,j; atom[n].tag=1; for (i=0; i<atom[n].numb; i++) { j=atom[n].bnd[i]; if (atom[j].tag==0) { atom[j].tag=1; tagatoms(atom,j); } } return; } /* dot product */ float dot(float *a,float *b) { return (a[0]*b[0]+a[1]*b[1]+a[2]*b[2]); } /* cross product */ void cross(float *a,float *b,float *c) { c[0]=a[1]*b[2]-a[2]*b[1]; c[1]=a[2]*b[0]-a[0]*b[2]; c[2]=a[0]*b[1]-a[1]*b[0]; return; } /* vector subtract */ void sub(float *a,float *b,float *c) { c[0]=a[0]-b[0]; c[1]=a[1]-b[1]; c[2]=a[2]-b[2]; return; } /* vector length */ float len(float *a) { return (sqrt(SQR(a[0])+SQR(a[1])+SQR(a[2]))); } /* return a random float between lo and hi */ float frand(float lo, float hi) { return (((float)random() / (float)MAXINT) * (hi - lo) + lo); } /* calculates a "chi" for fitting (not used yet ...) */ -(float)calcChi:sender { int i; int lb; float l,ch; float *a1,*a2; lb=mol[0].numlb; ch=0.0; for (i=0; i<lb; i++) { a1=(float *)&mol[0].coord[mol[0].lb[i].n1]; a2=(float *)&mol[0].coord[mol[0].lb[i].n2]; l=mol[0].lb[i].d; ch+=fabs(l-sqrt(SQR(a1[0]-a2[0])+SQR(a1[1]-a2[1])+SQR(a1[2]-a2[2]))); } ch/=lb; return ch; } /* fitting stuff ... */ -minStep:sender { int i,j,k,m,c,c2; float x,y,z,d,d2,v1,v2; float x1,y1,z1; float *a0,*a1; short *ba; c2=0; if (stepmode==0) { for (i=0; i<NSTEP0; i++) { x1=y1=z1=0.0; j=random()%mol[0].numa; k=mol[0].atom[j].numb; a0=(float *)&mol[0].coord[j]; for (m=0; m<k; m++) { a1=(float *)&mol[0].coord[mol[0].atom[j].bnd[m]]; x=a1[0]-a0[0]; y=a1[1]-a0[1]; z=a1[2]-a0[2]; d2=sqrt(SQR(x)+SQR(y)+SQR(z)); /* d=(1.0/mol[0].atom[j].bndl[m]-1.0/d2); */ d=1.0; x1+=x*d; y1+=y*d; z1+=z*d; } x1/=(float)k; y1/=(float)k; z1/=(float)k; a0[0]+=x1; a0[1]+=y1; a0[2]+=z1; } chi=[self calcChi:self]; c2=100; } else if (stepmode==1) { for (i=0; i<NSTEP; i++) { j=random()%mol[0].numa; m=1; c=0; ba=mol[0].atom[j].bnd; a0=(float *)&mol[0].coord[j]; v1=0.0; for (k=0; k<mol[0].numa; k++) { if (k==j||k==ba[0]||k==ba[1]||k==ba[2]||k==ba[3]) continue; a1=(float *)&mol[0].coord[k]; d=sqrt(SQR(a0[0]-a1[0])+SQR(a0[1]-a1[1])+SQR(a0[2]-a1[2])); v1+=1.0/d; } while (m&&c++<30) { x1=frand(-cstep,cstep); y1=frand(-cstep,cstep); z1=frand(-cstep,cstep); v2=0.0; a0[0]+=x1; a0[1]+=y1; a0[2]+=z1; m=0; for (k=0; k<mol[0].numa; k++) { if (k==j||k==ba[0]||k==ba[1]||k==ba[2]||k==ba[3]) continue; a1=(float *)&mol[0].coord[k]; d=sqrt(SQR(a0[0]-a1[0])+SQR(a0[1]-a1[1])+SQR(a0[2]-a1[2])); v2+=1.0/d; } /* if (v1<v2) m=1; else m=0; if (m) { a0[0]-=x1; a0[1]-=y1; a0[2]-=z1; continue; } */ d=v1/(mol[0].numa*6.0); v2/=mol[0].numa*6.0; x=[self calcChi:self]; if (x+v2>chi+d) { m=1; a0[0]-=x1; a0[1]-=y1; a0[2]-=z1; continue; } chi=x; } if (c<=30) c2++; } [chiDevDis setFloatValue:chi]; } [c2Dis setIntValue:c2]; if (c2<3) { cstep*=.9; [csDis setFloatValue:cstep]; } return self; } /* more fitting stuff */ -resetCs:sender { cstep=1.0; [csDis setFloatValue:cstep]; return self; } /* this too ... */ -togFit:sender { stepmode=[sender intValue]; return self; } -error:(char *)msg :(BOOL)xit { [errorMsg setStringValue:(void *)msg]; [errorPanel makeKeyAndOrderFront:self]; [errorPanel flushWindow]; if (xit) { sleep(5); exit(1); } return self; } @end
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.