This is NrgObj.m in view mode; [Download] [Up]
/* 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
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.