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;
}
@endThese are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.