/* NrgObj.m - Copyright 1993 Steve Ludtke */
/* This object is responsible for Amber energy calculations */
/* Sorry, I haven't gotten around to documenting it yet ... */

#include <ctype.h>

#import "Mtypes.h"
#import "NrgObj.h"
#import "MolObj.h"
#define A0 0.5292

float hydrogenic(int Z,float r);

@implementation NrgObj
-init
{
[super init];
bondl=NULL;
nonb=NULL;
hbond=NULL;
nmol=0;
rmax=100.0;
flag=NRG_HBOND|NRG_NONB|NRG_ELST;
return self;
}

-startup:(Molecule *)Mol :(struct ELINFO *)Elinfo
{
mol=Mol;
elinfo=Elinfo;
return self;
}

-readPotDat
{
char *bpath;
FILE *in;
char s[120];
unsigned char an,tp;
int i,j,k,l,n;

bpath=(char *)[[NXBundle mainBundle] directory];

sprintf(s,"%s/amb.nonbond.dat",bpath);
if ((in=fopen(s,"r"))==NULL) return self;
fgets(s,80,in);
sscanf(s," %d",&n);
nonb=malloc(sizeof(struct HBOND)*n);
for (i=0; i<n; i++) {
	fgets(s,80,in);
	j=[molobj parseAtom:s :&nonb[i].an :&nonb[i].tp];
	if (nonb[i].an==NEL) printf("Error in nonb line %d\n",i+2);
	sscanf(&s[j]," %f %f",&nonb[i].r,&nonb[i].e);
}
nnonb=n;
fclose(in);

sprintf(s,"%s/amb.hbond.dat",bpath);
if ((in=fopen(s,"r"))==NULL) return self;
fgets(s,80,in);
sscanf(s," %d",&n);
hbond=malloc(sizeof(struct HBOND)*n);
for (i=0; i<n; i++) {
	fgets(s,80,in);
	j=[molobj parseAtom:s :&an :&tp];
	if (an==NEL) printf("Error in bhond line %d 1\n",i+2);
	for (l=0; l<nnonb; l++) if (an==nonb[l].an&&tp==nonb[l].tp) break;
	if (l==nnonb) printf("Parse error %s\n",s);
	hbond[i].at1=l;
	k=[molobj parseAtom:&s[j] :&an :&tp];
	if (an==NEL) printf("Error in bhond line %d 2\n",i+2);
	for (l=0; l<nnonb; l++) if (an==nonb[l].an&&tp==nonb[l].tp) break;
	if (l==nnonb) printf("Parse error %s\n",s);
	hbond[i].at2=l;
	sscanf(&s[j+k]," %f %f",&hbond[i].c,&hbond[i].d);
}
nhbond=n;
fclose(in);

sprintf(s,"%s/amb.bondl.dat",bpath);
if ((in=fopen(s,"r"))==NULL) return self;
fgets(s,80,in);
sscanf(s," %d",&n);
bondl=malloc(sizeof(struct HBOND)*n);
for (i=0; i<n; i++) {
	fgets(s,80,in);
	j=[molobj parseAtom:s :&an :&tp];
	if (an==NEL) printf("Error in bondl line %d 1\n",i+2);
	for (l=0; l<nnonb; l++) if (an==nonb[l].an&&tp==nonb[l].tp) break;
	if (l==nnonb) printf("Parse error bondl 1 %d %s\n",i+2,s);
	bondl[i].at1=l;
	k=[molobj parseAtom:&s[j] :&an :&tp];
	if (an==NEL) printf("Error in bondl line %d 2\n",i+2);
	for (l=0; l<nnonb; l++) if (an==nonb[l].an&&tp==nonb[l].tp) break;
	if (l==nnonb) printf("Parse error bondl 1 %d %s\n",i+2,s);
	bondl[i].at2=l;
	sscanf(&s[j+k]," %f %f",&bondl[i].k,&bondl[i].r);
}
nbondl=n;
fclose(in);

if (nmol!=0) [self asnAmber];
return self;
}

/* assign amber atom and bond types */
-asnAmber
{
int i,j;
char at1,at2;

for (i=0; i<mol[0].numa; i++) {
	for (j=0; j<nnonb; j++) 
		if (mol[0].atom[i].anum==nonb[j].an && mol[0].atom[i].type==nonb[j].tp) break;
	if (j==nnonb) mol[0].atom[i].ambtp=-1;
	else mol[0].atom[i].ambtp=j;
}

for (i=0; i<mol[0].numlb; i++) {
	j=mol[0].lb[i].n1;
	at1=mol[0].atom[j].ambtp;
	j=mol[0].lb[i].n2;
	at2=mol[0].atom[j].ambtp;

	for (j=0; j<nbondl; j++) {
		if (at1==bondl[j].at1 && at2==bondl[j].at2) break;
		if (at2==bondl[j].at1 && at1==bondl[j].at2) break;
	}
	if (j==nbondl) mol[0].lb[i].ambtp=-1;
	else mol[0].lb[i].ambtp=j;
}

return self;
}

-(float)calcE:sender
{
int i,j,k,t1,t2;
float E=0.0,r,e,r0,t;

if (nmol==0) return 0.0;
if (bondl==NULL) [self readPotDat];

/* nonbonded */
for (i=0; i<mol[0].numa; i++) {
	for (j=i+1; j<mol[0].numa; j++) {
		for (k=0; k<4; k++) if (j==mol[0].atom[i].bnd[k]) break;
		if (k!=4) continue;
		r=sqrt(SQR(mol[0].coord[i][0]-mol[0].coord[j][0])+ 
			SQR(mol[0].coord[i][1]-mol[0].coord[j][1])+ 
			SQR(mol[0].coord[i][2]-mol[0].coord[j][2]));
		if (r>rmax) continue;
		if (r<.7) r=.7;
		t1=mol[0].atom[i].ambtp;
		t2=mol[0].atom[j].ambtp;

/* see if this should use HBOND potential or NONB potenetial */
		k=nhbond;
		if ((flag&NRG_HBOND)&&(t1<12&&t2<12)&&(t1<5 || t2<5)) {
			for (k=0; k<nhbond; k++) {
				if (t1==hbond[k].at1&&t2==hbond[k].at2) break;
				if (t2==hbond[k].at1&&t1==hbond[k].at2) break;
			}
		}

/* HBOND */
		if (k<nhbond) {
			E+=hbond[k].c/pow(r,12.0)-hbond[k].d/pow(r,10.0);
		/*	E+=t;
			printf("hb %d %d  %f %f\n",i,j,r,t); */
		}
/* NONB */
		else if (flag&NRG_NONB) {
			r0=(nonb[t1].r+nonb[t2].r)/2.0;
			e=sqrt(nonb[t1].e*nonb[t2].e);
			t=(e*(pow(r0/r,12.0)-2.0*pow(r0/r,6.0)));
			E+=t;
/*			if (t>30.0) printf("nb %d %d  %f %f   %f %f\n",i,j,r0,e,r,t);*/
		}
/* electrostatic */
		if (flag&NRG_ELST) {
			E+=332.0*mol[0].atom[i].charge*mol[0].atom[j].charge/SQR(r);
		/*	E+=t;
			printf("es %d %d  %f %f\n",i,j,r,t); */
		}
	}
}
		
return E;
}

-setFlag:sender
{
flag=0;
if ([[sender cellAt:0 :0] intValue]) flag|=NRG_HBOND;
if ([[sender cellAt:1 :0] intValue]) flag|=NRG_NONB;
if ([[sender cellAt:2 :0] intValue]) flag|=NRG_ELST;
return self;
}

-setRmax:sender
{
rmax=[sender floatValue];
if (rmax<1.0) { rmax=1.0; [sender setFloatValue:rmax]; }
return self;
}

-recalc:sender
{
[Edisp setFloatValue:[self calcE:self]];
return self;
}

-update:(int)Nmol level:(int)level
{
if (nmol!=Nmol || level>=U_BONDS) {
	nmol=Nmol;
	[self asnAmber];
}
return self;
}	

-series:sender
{
float x,y,E;
FILE *out;
char s[120];

sprintf(s,"%s/prot.nrg",getenv("HOME"));
out=fopen(s,"w");
for (y=-180.0; y<=180.0; y+=10.0) {
	for (x=-180.0; x<=180.0; x+=10.0) {
		[molobj setPsi:x];
		[molobj setPhi:y];
		E=[self calcE:self];
		fprintf(out,"%f\t%f\t%f\n",x,y,E);
		[Edisp setFloatValue:E];
		[[scanstat cellAt:1 :0] setFloatValue:x];
		[[scanstat cellAt:2 :0] setFloatValue:y];
	}
}
fclose(out);

return self;
}

-series2:sender
{
float x,y,z,E,e2,z2;
FILE *out,*out2;
char s[120];

sprintf(s,"%s/prot.nrg",getenv("HOME"));
out=fopen(s,"w");
sprintf(s,"%s/prot.chi",getenv("HOME"));
out2=fopen(s,"w");
for (y=-70.0; y<=20.0; y+=6.0) {
	for (x=-110.0; x<=-20.0; x+=6.0) {
		[molobj setPsi:x];
		[molobj setPhi:y];
		e2=MAXFLOAT;
		z2=-180.0;
		for (z=-180.0; z<180.0; z+=20.0) {
			[molobj setChi:z];
			E=[self calcE:self];
			if (E<e2) { e2=E; z2=z; } 
			[[scanstat cellAt:3 :0] setFloatValue:z];
		}
		E=e2;
		fprintf(out,"%f\t%f\t%f\n",x,y,E);
		fprintf(out2,"%f\t%f\t%f\n",x,y,z2);
		[Edisp setFloatValue:E];
		[[scanstat cellAt:1 :0] setFloatValue:x];
		[[scanstat cellAt:2 :0] setFloatValue:y];
	}
}
fclose(out);
fclose(out2);

return self;
}

-calcElecDenProf:sender
{
float x,y,z,d,e;
FILE *out;

out=fopen("/users/steve/x.out","w");
if (out==NULL) return self;
for (y=-20.0; y<=20.0; y+=.05) {
	d=0.0;
/*
	for (z=-10.0; z<=10.0; z+=.5) {
		for (x=-10.0; x<=10.0; x+=.5) {
			e=[self elecDen:x :y :z];
			if (e>50.0) e=[self elecDen:x+.05 :y+.05 :z+.05];
			if (e>50.0) e=[self elecDen:x :y :z];
			d+=e;
		}
	}
*/
	d=[self elecDenY:y];
	fprintf(out,"%f\t%f\n",y,d);
}
fclose(out);
return self;
}

-(float)elecDenY:(float)y
{
int i;
float d=0.0,r;

for (i=0; i<mol[0].numa; i++) {
	r=fabs(y-mol[0].coord[i][1]);
	d+=hydrogenic(mol[0].atom[i].anum+1,r);
}

return d;
}


-(float)elecDen:(float)x :(float)y :(float)z
{
int i;
float d=0.0,r;

for (i=0; i<mol[0].numa; i++) {
	r=sqrt(SQR(x-mol[0].coord[i][0])+SQR(y-mol[0].coord[i][1])+SQR(z-mol[0].coord[i][2]));
	d+=hydrogenic(mol[0].atom[i].anum+1,r);
}

return d;
}

/* approximates the electron density of an atom at a radius R */
/* uses only the radial wavefunctions (s orbitals) */
float hydrogenic(int z,float r)
{
float d=0.0,zap,sig,t,Z;
int i;

Z=z;
zap=pow(Z/A0,3.0)/M_PI;
sig=Z*r/A0;

/* K */
t=zap*exp(-sig*2.0);
for (i=0; i<2&&i<z; i++) {
	d+=t;
}

/* L */
t=zap/32.0*pow(2.0-sig,2.0)*exp(-sig);
for (i=2; i<10&&i<z; i++) {
	d+=t;
}

/* M */
t=zap/19683.0*pow(27-18*sig+2*sig*sig,2.0)*exp(-sig*.6666667);
for (i=10; i<54&&i<z; i++) {
	d+=t;
}
return(d);
}
@end
