/* 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