steric_1.11
|
Makefile,
Makefile.sgi,
README.steric,
contplot,
craig.c,
craig.h,
crystal.c,
crystal.h,
integrat.c,
integrat.h,
leach.c,
leach.h,
long_steric,
makeit,
mapcont,
mapcont.c,
mapprof,
mapprof.c,
profplot,
proja.c,
proja.h,
ryan.c,
ryan.h,
ryan_perm.c,
ryan_perm.h,
ryan_quad.c,
ryan_quad.h,
steraid.c,
steraid.h,
stercalc.c,
stercalc.h,
stercomm.c,
stercomm.h,
sterdefn.h,
stererr.h,
sterfile.c,
sterfile.h,
stergrap.c,
stergrap.h,
stergrp.c,
stergrp.h,
steric,
steric.TeX,
steric.err,
steric.grp,
steric.hlp,
steric.ini,
steric.par,
stermain.c,
stermem.c,
stermem.h,
sterover.c,
sterover.h,
sterplot,
stertext.c,
stertext.h,
test.bgf,
test.inp,
vectors.c,
vectors.h,
|
|
|
/**************************************************************************/
/****************** Craig solid angle calculations **********************/
/**************************************************************************/
#include
#include
#include
#include
#include "sterdefn.h"
#include "integrat.h"
#include "stermem.h"
#include "steraid.h"
#include "stercalc.h"
#include "stertext.h"
#include "craig.h"
/**************************************************************************/
/****************** ellipse memory management ***************************/
/**************************************************************************/
void Initialize_Ellipse(Elps *elps, Atms *at)
{
elps->a=0.0; /* long axis */
elps->b=0.0; /* short axis */
elps->c=0.0; /* shift off origin */
elps->delta=0.0; /* atom angle from vector G */
elps->A=at; /* pointer to atom memory */
elps->I[0]=Vequal(0.0,0.0,0.0);
elps->I[1]=Vequal(0.0,0.0,0.0);
elps->phi[0]=PI/2;
elps->phi[1]=PI/2; /* angles of intersepts from main axis */
elps->mode=0; /* for special ellipse modes */
}
/**************************************************************************/
void Initialize_Over(Over *over, Atms *A[MAX_OVER], int order)
{
int i=0;
if(order>MAX_OVER) order=MAX_OVER;
if(over!=NULL)
{
for(i=0;iEl[i],A[i]);
over->G=Vequal(0.0,0.0,0.0); /* vector to center of overlap */
over->order=order; /* overlap order (eg. 3=triple */
over->nell=order; /* number of ellipses bounding over */
over->solid=0.0; /* solid angle calculated */
}
/* next and prev pointers have been set so don't touch */
}
/**************************************************************************/
void Copy_Ellipse(Elps *new, Elps *old)
{
new->a=old->a;
new->b=old->b;
new->c=old->c;
new->delta=old->delta;
new->I[0]=VequalV(old->I[0]);
new->I[1]=VequalV(old->I[1]);
new->phi[0]=old->phi[0];
new->phi[1]=old->phi[1];
}
/**************************************************************************/
void Copy_Over(Over *new, Over *old)
{
int i=0;
if((new!=NULL)&&(old!=NULL))
{
for(i=0;iEl[i],&old->El[i]);
new->G=VequalV(old->G);
new->order=old->order;
new->nell=old->nell;
new->solid=old->solid;
}
}
/**************************************************************************/
/**************************************************************************/
Over *New_Over(Over *old)
{
Over *new=NULL;
if((new=(Over *)malloc(sizeof(Over)))==NULL)
{
Error_Message(E_NOMEM,"New Over");
return(NULL);
}
new->next=NULL;
new->prev=NULL;
if(old==NULL) return(new);
new->next=old->next;
new->prev=old;
if(old->next) old->next->prev=new;
old->next=new;
return(new);
}
/**************************************************************************/
/**************************************************************************/
Over *First_Over(Over *current)
{
if(current==NULL) return(NULL);
while(current->prev!=NULL) current=current->prev;
return(current);
}
/**************************************************************************/
Over *Last_Over(Over *current)
{
if(current==NULL) return(NULL);
while(current->next!=NULL) current=current->next;
return(current);
}
/**************************************************************************/
/**************************************************************************/
int Get_Over_Number(Over *current)
{
int count=0;
if(current==NULL) return(0);
while(count++,current->prev!=NULL) current=current->prev;
return(count);
}
/**************************************************************************/
/**************************************************************************/
Over *Goto_Over(Over *current, int num)
{
int count=0;
current=First_Over(current);
if(current==NULL) return(NULL);
while(count++,(current->next!=NULL)&&(count!=num)) current=current->next;
return(current);
}
/**************************************************************************/
/**************************************************************************/
Over *Next_Over(Over *current)
{
if(current==NULL) return(NULL);
if(current->next==NULL) return(current);
return(current->next);
}
/**************************************************************************/
Over *Previous_Over(Over *current)
{
if(current==NULL) return(NULL);
if(current->prev==NULL) return(current);
return(current->prev);
}
/**************************************************************************/
/**************************************************************************/
Over *Close_Over(Over *current)
{
Over *old;
if(current==NULL) return(NULL);
if(current->prev==NULL)
{
if(current->next==NULL)
{
free(current);
return(NULL);
}
current=current->next;
free(current->prev);
current->prev=NULL;
return(current);
}
if(current->next==NULL)
{
current=current->prev;
free(current->next);
current->next=NULL;
return(current);
}
old=current;
current=old->prev;
old->prev->next=old->next;
old->next->prev=old->prev;
free(old);
return(current);
}
/**************************************************************************/
void Close_All_Overlaps(Over *over, char *title, unsigned short mode)
{
Over *o=NULL;
char line[LINELEN];
int count=0;
for(o=First_Over(over),count=0;o!=NULL;count++,o=o->next);
if((count>0)&&(mode&VIS_BIT))
{
sprintf(line,"Closing temporary overlap calculation memory for %d %s overlaps",count,title);
Out_Message(line,O_NEWLN);
}
for(over=First_Over(over);over!=NULL;over=Close_Over(over));
}
/**************************************************************************/
/****************** ellipse calculations ********************************/
/**************************************************************************/
void Find_2atom_Intersepts(double alpha, double beta
,Vector V[2], Vector I[2])
{
int i=0;
double D=0.0,T=1.0;
double p=0.0,q=0.0,l=0.0,m=0.0;
double a=0.0,b=0.0,c=0.0;
while((i<3)&&(AlmostZero(T= V[0].x * V[1].y - V[1].x * V[0].y)))
{
Rotate_Indices_Right(&V[0]);
Rotate_Indices_Right(&V[1]);
i++;
}
p= ( V[0].x * cos(beta) - V[1].x * cos(alpha))/T;
q= ( V[1].x * V[0].z - V[0].x * V[1].z)/T;
l= cos(alpha) / V[0].x - p * V[0].y / V[0].x ;
m= V[0].z / V[0].x + q * V[0].y / V[0].x ;
a= m*m + q*q + 1 ;
b= 2*p*q - 2*l*m ;
c= l*l + p*p - 1 ;
D= b*b -4*a*c;
if(D<0.0)
{
I[0]=Vequal(0.0,0.0,0.0);
I[1]=Vequal(0.0,0.0,0.0);
}
else
{
I[0].z = (-b-sqrt(D))/(2*a);
I[0].x = l-m*I[0].z;
I[0].y = p+q*I[0].z;
I[1].z = (-b+sqrt(D))/(2*a);
I[1].x = l-m*I[1].z;
I[1].y = p+q*I[1].z;
I[0]=unit_vector(I[0]);
I[1]=unit_vector(I[1]);
}
for(;i>0;i--)
{
Rotate_Indices_Left(&V[0]);
Rotate_Indices_Left(&V[1]);
Rotate_Indices_Left(&I[0]);
Rotate_Indices_Left(&I[1]);
}
}
/**************************************************************************/
void Find_2atom_G(Over *over, Vector V[2])
{
double gap=0.0;
Find_2atom_Intersepts(over->El[0].A->SVangle,over->El[1].A->SVangle
,V,over->El[0].I);
over->El[1].I[0]=VequalV(over->El[0].I[0]);
over->El[1].I[1]=VequalV(over->El[0].I[1]);
over->G=unit_vector(Vsum(over->El[0].I[0],over->El[0].I[1]));
over->El[0].delta=VangleV(over->G,over->El[0].A->v);
over->El[1].delta=VangleV(over->G,over->El[1].A->v);
gap=over->El[0].delta+over->El[1].delta
-VangleV(over->El[0].A->v,over->El[1].A->v);
if((gap>0.0)&&!(Small(gap))) /* fix to solve inadequacies in ellipse */
{ /* parameter calculations */
if(over->El[0].deltaEl[1].delta) over->El[0].mode|=S_NEGS;
else over->El[1].mode|=S_NEGS;
}
}
/**************************************************************************/
void Get_Ellipse_Mode(Elps *Elly, Vector G)
{
Vector A={0.0,0.0,0.0}; /* vector to intersept average */
double angle=0.0; /* angle between G and A */
/* work out whether ellipse is +ve or -ve */
A=SxV(0.5,Vsum(Elly->I[0],Elly->I[1]));
angle=VangleV(A,G);
if(!AlmostZero(angle))
{
if(VangleV(Elly->A->v,A)mode|=S_NEGS;
}
}
/**************************************************************************/
void Calculate_Ellipse_Intersept_Angles(Elps *Elly, double theta, double phi)
{
double angle=0.0;
Vector X[2]; /* rotated intersept vectors */
Vector A={0.0,0.0,0.0}; /* vector to atom rotated */
int i=0;
/* calculate intersept vectors in plane normal to vector G */
A=Rotate_About_Z(Elly->A->v,-1.0*phi);
A=Rotate_About_Y(A,-1.0*theta);
A.z=0.0;
for(i=0;i<2;i++)
{
X[i]=Rotate_About_Z(Elly->I[i],-1.0*phi);
X[i]=Rotate_About_Y(X[i],-1.0*theta);
X[i].z=0.0;
}
/* calculate intersept angles phi[] */
angle=VangleV(X[0],X[1]);
Elly->phi[0]=VangleV(X[0],A);
Elly->phi[1]=VangleV(X[1],A);
if(AlmostZero(2*PI-angle-Elly->phi[0]-Elly->phi[1]))
{ /* normal situation with +ve and -ve phi's */
Elly->phi[0]=PI-Elly->phi[0];
Elly->phi[1]=PI-Elly->phi[1];
}
else
{ /* both phi's on same side */
if(Elly->phi[0]phi[1])
{
Elly->phi[0]=PI-Elly->phi[0];
Elly->phi[1]-=PI;
}
else
{
Elly->phi[1]=PI-Elly->phi[1];
Elly->phi[0]-=PI;
}
}
}
/**************************************************************************/
void Calculate_Ellipse_Parameters(Elps *Elly)
{
double s=0.0;
/* calculate ellipse coefficients using the Leach equations (White et al.)*/
s=cos(2*Elly->A->SVangle)+cos(2*Elly->delta);
Elly->a=sin(2*Elly->A->SVangle)/s;
Elly->b=msqrt(2*pow(sin(Elly->A->SVangle),2)/s,"msqrt -> Calculate Ellipse Parameters (b)");
Elly->c=sin(2*Elly->delta)/s;
}
/**************************************************************************/
void Setup_Two_Ellipses(Over *over, Atms *A[MAX_OVER], char mode)
{
Vector V[2]={{0.0,0.0,0.0},{0.0,0.0,0.0}};
double theta=0.0,phi=0.0;
switch(mode)
{
case FIRST : over->G=unit_vector(A[0]->v);
over->nell=1;
over->El[0].A=A[0];
break;
case SECOND : over->G=unit_vector(A[1]->v);
over->nell=1;
over->El[0].A=A[1];
break;
case BOTH : V[0]=unit_vector(over->El[0].A->v);
V[1]=unit_vector(over->El[1].A->v);
Find_2atom_G(over,V);
theta=Get_Vector_Theta(over->G);
phi=Get_Vector_Phi(over->G);
over->El[0].delta=VangleV(over->El[0].A->v,over->G);
over->El[1].delta=VangleV(over->El[1].A->v,over->G);
Calculate_Ellipse_Intersept_Angles(&over->El[0],theta,phi);
Calculate_Ellipse_Intersept_Angles(&over->El[1],theta,phi);
Calculate_Ellipse_Parameters(&over->El[0]);
Calculate_Ellipse_Parameters(&over->El[1]);
break;
default : break;
}
}
/**************************************************************************/
/****************** integration calculations ****************************/
/**************************************************************************/
static double iA;
static double iB;
static double iC;
static int isign;
void Prepare_Ellipse(Elps *Elly)
{
iA=Elly->a;
iB=Elly->b;
iC=Elly->c;
if(Elly->mode&S_NEGS)
{
isign=-1;
printf("n");
}
else isign=1;
}
/**************************************************************************/
double Ellipse(double phi) /* this is the general ellipse equation */
{
double p,q;
double s;
double eta;
p=cos(phi)*iC/(iA*iA);
q=pow(cos(phi)/iA,2)+pow(sin(phi)/iB,2);
s=p*p+q*(1-(iC/iA)*(iC/iA));
if (s<0.0) s=0.0;
eta=(msqrt(s,"msqrt -> Ellipse (eta)")-isign*p)/q;
return(1/msqrt(1+eta*eta,"msqrt -> Ellipse"));
}
/**************************************************************************/
/**************************************************************************/
/****************** solid angle calculations ****************************/
/**************************************************************************/
/**************************************************************************/
void Add_Phi(Over *over, int order, double excess)
{
int count,n;
double phi=0.0;
while(excess>0.0)
{
for(count=0,n=0;countEl[count].phi[1]El[count].phi[1]>phi) phi=over->El[count].phi[1];
}
if(over->El[count].phi[0]El[count].phi[0]>phi) phi=over->El[count].phi[0];
}
}
phi=PI/2.0-phi;
if(excess/(double)n<=phi) phi=excess/(double)n;
for(count=0;countEl[count].phi[0]>0.0) over->El[count].phi[0]-=phi;
if(over->El[count].phi[1]>0.0) over->El[count].phi[1]-=phi;
}
excess-=phi*(double)n;
if(AlmostZero(excess)) break;
}
}
/**************************************************************************/
void Cut_Phi(Over *over, int order, double excess)
{
int count,n;
double phi=100.0;
while(excess>0.0)
{
for(count=0,n=0;countEl[count].phi[1]>0.0)
{
n++;
if(over->El[count].phi[1]El[count].phi[1];
}
if(over->El[count].phi[0]>0.0)
{
n++;
if(over->El[count].phi[0]El[count].phi[0];
}
}
if(excess/(double)n<=phi) phi=excess/(double)n;
for(count=0;countEl[count].phi[0]>0.0) over->El[count].phi[0]-=phi;
if(over->El[count].phi[1]>0.0) over->El[count].phi[1]-=phi;
}
excess-=phi*(double)n;
if(AlmostZero(excess)) break;
}
}
/**************************************************************************/
/****************** single atom calculations ****************************/
/**************************************************************************/
double Get_All_Single_Solids(Mol *M, unsigned mode)
{
int count=0;
double solid=0.0;
Atms *A=NULL;
char line[LINELEN];
for(A=M->atoms,count=0;A!=NULL;A=A->next)
{
count++; if(count>M->num_atoms) Error_Message(E_ATMNUM,"Get All Single Solids");
if(A->stat&CALC_BIT)
{
if(mode&VIS_BIT) Out_Message("+ ",O_BLANK);
solid+=Single_Atom_Solid_Angle(A,mode);
if(mode&VIS_BIT)
{
sprintf(line,"= %f",solid);
Out_Message(line,O_NEWLN);
}
A->count=1;
}
else A->count=0;
}
return(solid);
}
/**************************************************************************/
/**************** two atom overlap calculation **************************/
/**************************************************************************/
double Double_Overlap_Solid(Over *over, double eps, unsigned mode)
{
double totphi=0.0;
double integral=0.0;
int count=0;
char line[150];
/* check for excessively large phi range */
for(count=0;count<2;count++)
totphi+=over->El[count].phi[0]+over->El[count].phi[1];
totphi-=2.0*PI;
if(!AlmostZero(totphi))
{
if(mode&VIS_BIT)
{
sprintf(line,"Double Overlap Solid [%7.4f]",totphi);
Error_Message(E_PNOTPI,line);
}
if(totphi>0.0) Cut_Phi(over,2,totphi);
if(totphi<0.0) Add_Phi(over,2,-1.0*totphi);
}
/* calculate right ellipse */
Prepare_Ellipse(&over->El[0]);
integral+=qsimp(Ellipse,0.0,over->El[0].phi[1],eps);
/* calculate left ellipse */
Prepare_Ellipse(&over->El[1]);
integral+=qsimp(Ellipse,0.0,over->El[1].phi[1],eps);
/* calculate solid angle */
integral=2*PI-2*integral;
/* display answers if in normal total calculation */
if(mode&VIS_BIT)
{
sprintf(line,"double(%s[%d]-%s[%d]):[%f] "
,over->El[0].A->name,Get_Atom_Number(over->El[0].A)
,over->El[1].A->name,Get_Atom_Number(over->El[1].A),integral);
Out_Message(line,O_BLANK);
}
if(integral<0.0)
{
Error_Message(E_NEGSLD,"Double Overlap Solid");
return(0.0);
}
return(integral);
}
/**************************************************************************/
double Remove_All_Double_Overlaps(Mol *M, Over **current, double eps
,unsigned mode, double solid)
{
int count,i;
Over *new=NULL;
Over *over=NULL;
double chi=0.0;
Atms *A[MAX_OVER];
char line[LINELEN];
for(i=0;iatoms,count=0;A[0]!=NULL;A[0]=A[0]->next)
{
count++; if(count>M->num_atoms) Error_Message(E_ATMNUM,"Remove All Double Overlaps");
if((A[0]->SVangle!=0.0)&&(A[0]->stat&CALC_BIT))
{ /* count A[1] though whole list */
for(A[1]=A[0]->next,i=0;A[1]!=NULL;A[1]=A[1]->next)
{
i++; if(i>M->num_atoms) Error_Message(E_ATMNUM,"Remove All Double Overlaps");
if ((A[0]!=A[1])&&(A[1]->SVangle!=0.0)&&(A[1]->stat&CALC_BIT))
{
chi=VangleV(A[0]->v,A[1]->v);
if (chi<(A[0]->SVangle+A[1]->SVangle))
{ /* A[0]-A[1] steric overlap found */
if((new=New_Over(over))!=NULL)
{
over=new; /* create the memory for counted doubles */
Initialize_Over(over,A,2);
if(mode&VIS_BIT) Out_Message("- ",O_BLANK);
if (chi<=fabs(A[0]->SVangle-A[1]->SVangle))
{ /* total overlap found */
Error_Message(E_UNSING,"Remove All Double Overlaps");
if(mode&VIS_BIT) Out_Message("- ",O_BLANK);
if (A[0]->SVangle>=A[1]->SVangle)
{ /* A[0] overlaps A[1] */
Setup_Two_Ellipses(over,A,SECOND);
solid-=Single_Atom_Solid_Angle(A[1],mode);
}
else
{ /* A[1] overlaps A[0] */
Setup_Two_Ellipses(over,A,FIRST);
solid-=Single_Atom_Solid_Angle(A[0],mode);
}
if(mode&VIS_BIT)
{
sprintf(line,"= %f",solid);
Out_Message(line,O_NEWLN);
}
}
else /* partial overlap found */
{
Setup_Two_Ellipses(over,A,BOTH);
solid-=Double_Overlap_Solid(over,eps,mode);
if(mode&VIS_BIT)
{
sprintf(line,"= %f",solid);
Out_Message(line,O_NEWLN);
}
}
}
}
}
}
}
}
*current=over;
return(solid);
}
/**************************************************************************/
/**************************************************************************/
/**************** multiple atom overlap calculation *********************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/** test for case where one atom is enveloped by another atom **/
int Inside_Atom(Atms *A, Atms *B)
{
if(A==B) return(1);
if(A->SVangle+VangleV(A->v,B->v)<=B->SVangle) return(1);
return(0);
}
/**************************************************************************/
/** test for that atom is not enveloped by any other atom **/
int Outside_All_Atoms(Atms *A)
{
Atms *atoms;
if(A==NULL) return(0);
for(atoms=First_Atom(A);atoms!=NULL;atoms=atoms->next)
{
if(atoms==A) atoms=atoms->next;
if(atoms==NULL) break;
if(Inside_Atom(A,atoms)) return(0);
}
return(1);
}
/**************************************************************************/
/** test for case where overlap region in entirely within an atom **/
int Atom_Covers_Overlap(Atms *A, Over *O)
{
int i;
for(i=0;inell;i++)
{
if(VangleV(O->El[i].I[0],A->v)>A->SVangle) return(0);
if(VangleV(O->El[i].I[1],A->v)>A->SVangle) return(0);
}
return(1);
}
/**************************************************************************/
/* find overlap region that doesn't involve specified atom **/
int Find_Over_Without_Atom(Over *po[MAX_OVER], Atms *at, int o[2]
,int order, int start)
{
int i,k,a;
for(i=start;iEl[k].A) a++;
}
if(a>1) return(0); /* found same atom twice ?? */
if(a==0)
{
if(at==NULL) return(0);
else
{
o[1]=i; /* found it ! */
return(1);
}
}
}
return(0);
}
/**************************************************************************/
/* swap positions of two pointers to overs in the array po */
void Swap_Overs(Over *po[MAX_OVER], int o[2])
{
Over *temp;
temp=po[o[1]];
po[o[1]]=po[o[0]];
po[o[0]]=temp;
}
/**************************************************************************/
/* order the n* "(n-1) overlaps" to facilitate discovery of nth overlap */
int Order_Previous_Overlaps(Over *po[MAX_OVER], Atms *at[MAX_OVER], int order)
{
Atms *atom=NULL;
int i,k,o[2],a=0;
if(order<3) return(0); /* only triple or greater considered */
for(i=0;iEl[i-1].A;
if(at[i]==NULL) return(0);
} /* find appropriate overlap for at[1] */
if(!Find_Over_Without_Atom(po,at[1],o,order,1)) return(0);
o[0]=1;
Swap_Overs(po,o); /* swap over 1 and o[1] */
for(k=0;kEl[k].A;
if(atom==NULL) return(0);
if(at[1]==atom) return(0); /* Find_Over_Without_Atom() failed */
for(i=2;i1) return(0); /* found same atom twice ?? */
if(a==0)
{
if(at[0]==NULL) at[0]=atom;
else return(0); /* two unknowns found */
}
}
/* order the rest of the overlaps */
for(i=2;iEl[k].A==at[a]) break;
}
if(a==order) return(0);
}
}
return(1);
}
/**************************************************************************/
/* find the correct atomic intersepts for particular pair of atoms */
int Find_Atomic_Intersepts(Over *dover, Atms *A[MAX_OVER]
,int a, int b
,Vector I[2], int inter[2])
{
int i;
double angle;
I[0]=Vequal(0.0,0.0,0.0);
I[1]=Vequal(0.0,0.0,0.0);
inter[0]=0;
inter[1]=0;
for(dover=First_Over(dover);dover!=NULL;dover=dover->next)
{ /* find ab overlap */
if(((dover->El[0].A==A[a])&&(dover->El[1].A==A[b]))
||((dover->El[1].A==A[a])&&(dover->El[0].A==A[b])))
{
inter[0]=1; /* check that both intersepts are within all atoms */
inter[1]=1;
for(i=0;((iEl[0].I[0],A[i]->v)-A[i]->SVangle;
if((angle>0.0)&&(!AlmostZero(angle))) inter[0]=0;
angle=VangleV(dover->El[0].I[1],A[i]->v)-A[i]->SVangle;
if((angle>0.0)&&(!AlmostZero(angle))) inter[1]=0;
}
if(inter[0]) I[0]=VequalV(dover->El[0].I[0]);
if(inter[1]) I[1]=VequalV(dover->El[0].I[1]);
return(1);
}
}
return(0);
}
/**************************************************************************/
/* determine the nature and parameters of the new nth order overlap */
Over *New_Multiple_Overlap(Over *over, Over *dover, int order
,Over *O[MAX_OVER], Atms *A[MAX_OVER])
{
int i=0,l=0,k=0,K=0,n=0,p=0;
Vector I[MAX_OVER-1][2]; /* vectors for phi[] calc. */
int inter[MAX_OVER-1][2];
char mode=S_NORMAL; /* normal overlap combination expected */
int multi_atom=0;
double theta=0.0,phi=0.0; /* G's angular position */
for(i=0;iorder=order;
return(over);
}
}
for(i=0;iEl[i].A=A[i]; /* set ellipse atom */
for(l=0;l=order) break;
Find_Atomic_Intersepts(dover,A,i,k,I[n],inter[n]);
if(inter[n][0]) p++;
if(inter[n][1]) p++;
}
if(p<2) return(Close_Over(over)); /* no multiple overlap found */
if((p!=2)&&(p!=2*n)) return(Close_Over(over)); /* unknown combo */
else if(p==2)
{ /* normal type of intersept */
for(K=0,n=0;KEl[i].I[n]=VequalV(I[K][0]);
n++;
}
if(n>2) return(Close_Over(over)); /* unexpected error */
if(inter[K][1])
{
over->El[i].I[n]=VequalV(I[K][1]);
n++;
}
if(n>2) return(Close_Over(over)); /* unexpected error */
}
if(n!=2) return(Close_Over(over)); /* unexpected error */
}
else if(K==2*n)
{ /* special type of intersept */
mode=S_GCENT;
multi_atom=i; /* atom i used for partial single atom */
over->El[i].mode|=S_GCENT;
}
}
/* calculate the vector G within the n-overlap region */
if(mode&S_GCENT)
{
over->G=VequalV(unit_vector(over->El[multi_atom].A->v));
}
else
{
over->G=Vequal(0.0,0.0,0.0);
for(i=0;iG=Vsum(over->G,over->El[i].I[0]);
over->G=Vsum(over->G,over->El[i].I[1]);
}
over->G=unit_vector(over->G);
}
/* check that G is not too near any atom centres */
for(i=0;iEl[i].delta=VangleV(over->G,over->El[i].A->v);
if((!(mode&S_GCENT))
&&(Similar(over->El[i].A->SVangle
,over->El[i].A->SVangle-over->El[i].delta
,S_PERC)))
{
over->El[i].mode|=S_GCENT;
over->G=unit_vector(A[i]->v);
mode|=S_GCENT;
}
}
/* make absolutely sure that this vector is indeed within all atoms */
for(i=0;iEl[i].delta=VangleV(over->G,over->El[i].A->v);
if(over->El[i].delta>A[i]->SVangle)
{
Error_Message(E_BAD_G,"New Multiple Overlap");
return(Close_Over(over));
}
}
/* calculate the ellipse parameters for all necessary ellipses */
for(i=0;iEl[i],over->G);
if(!(over->El[i].mode&S_GCENT))
{
theta=Get_Vector_Theta(over->G);
phi=Get_Vector_Phi(over->G);
Calculate_Ellipse_Intersept_Angles(&over->El[i],theta,phi);
Calculate_Ellipse_Parameters(&over->El[i]);
}
}
if(mode&S_GCENT)
{ /* centred ellipses dealt with differently */
for(i=0;iEl[i].mode&S_GCENT)
{
over->El[i].phi[0]=PI;
over->El[i].phi[1]=PI;
for(n=0;nEl[i].phi[0]-=over->El[n].phi[0];
over->El[i].phi[1]-=over->El[n].phi[1];
}
}
}
}
}
return(over);
}
/**************************************************************************/
/* perform the necessary integration on the new nth order overlap */
double Multiple_Overlap_Solid(Over *over, int order
,double eps, unsigned mode)
{
double totphi=0.0;
double integral=0.0;
double multi_solid=0.0; /* extra solid angle for S_GCENT mode */
int count;
char line[150];
/* check for excessively large phi range */
for(count=0;countEl[count].phi[0]+over->El[count].phi[1];
totphi-=2.0*PI;
if(!AlmostZero(totphi))
{
if(mode&VIS_BIT)
{
sprintf(line,"Multiple Overlap Solid [%7.4f]",totphi);
Error_Message(E_PNOTPI,line);
}
if(totphi>0.0) Cut_Phi(over,order,totphi);
if(totphi<0.0) Add_Phi(over,order,-1.0*totphi);
for(count=0;countEl[count].phi[0]+over->El[count].phi[1];
totphi-=2.0*PI;
}
/* integrate over each of the ellipses present */
if(AlmostZero(totphi))
{
for(count=0;countEl[count].mode&S_GCENT)
{
printf("m");
totphi=over->El[count].phi[0]+over->El[count].phi[1];
multi_solid=totphi/(2*PI)
*Single_Atom_Solid_Angle(over->El[count].A,0);
}
else
{
Prepare_Ellipse(&over->El[count]);
if(over->El[count].phi[0]<0.0)
{
integral+=qsimp(Ellipse,-over->El[count].phi[0]
,over->El[count].phi[1],eps);
}
else if(over->El[count].phi[1]<0.0)
{
integral+=qsimp(Ellipse,-over->El[count].phi[1]
,over->El[count].phi[0],eps);
}
else
{
integral+=qsimp(Ellipse,0.0,over->El[count].phi[0],eps);
integral+=qsimp(Ellipse,0.0,over->El[count].phi[1],eps);
}
}
}
/* calculate solid angle from total integration results */
integral=2*PI-totphi-integral+multi_solid;
}
/* display answers if in normal total calculation */
if(mode&VIS_BIT)
{
sprintf(line,"[%f] ",integral);
Out_Message(line,O_BLANK);
}
if(integral<0.0)
{
Error_Message(E_NEGSLD,"Multiple Overlap Solid");
return(0.0);
}
return(integral);
}
/**************************************************************************/
/* find nth order overlap and calculate the relevant solid angle */
double Calc_Multiple_Overlap(Over *po[MAX_OVER], Over **over
,Over *dover, int order
,double eps, unsigned mode
,char *multiple, double solid)
{
Over *new=NULL;
char line[LINELEN];
char atmstr[25];
Atms *A[MAX_OVER];
Over *O[MAX_OVER];
int i;
short sign=1;
/* initialize */
if(order/2*2==order) sign=1; else sign=-1;
for(i=0;i0) sprintf(line,"+ %s[%d](",multiple,new->nell);
else sprintf(line,"- %s[%d](",multiple,new->nell);
for(i=0;(iEl[i].A!=NULL);i++)
{
sprintf(atmstr,"%s[%d]-",A[i]->name,Get_Atom_Number(A[i]));
strcat(line,atmstr);
}
line[strlen(line)-1]=0;
strcat(line,") ");
Out_Message(line,O_BLANK);
}
if(new->solid==0.0)
solid+=sign*Multiple_Overlap_Solid(new,new->nell,eps,mode);
else solid+=sign*new->solid;
if(mode&VIS_BIT)
{
sprintf(line,"= %f",solid);
Out_Message(line,O_NEWLN);
}
}
}
return(solid);
}
/**************************************************************************/
/******* recursive routine to search possible (n-1) combinations **********/
/**************************************************************************/
double Multiple_Overlap(Over *pover, Over **over, Over *dover
,int order, int level, double eps
,unsigned mode, char *multiple, double solid)
{
static Over *po[MAX_OVER];
int i;
if(pover==NULL) return(solid); /* miss call */
if(level==0) /* initialize */
{
while(pover->prev!=NULL) pover=pover->prev;
for(i=0;iorder)
{
Error_Message(E_LEVHI,"Multiple Overlap");
return(solid); /* miss call */
}
if(pover->order!=order)
{
Error_Message(E_BADORD,"Multiple Overlap");
return(solid); /* miss call */
}
for(po[level]=pover;po[level]!=NULL;pover=pover->next,po[level]=pover)
{
if(levelnext,over,dover,order,1+level,eps,mode,multiple,solid);
}
else /* final depth, calculate */
{
solid=Calc_Multiple_Overlap(po,over,dover,order,eps,mode,multiple,solid);
}
}
return(solid);
}
/**************************************************************************/
/**************************************************************************/
/****** Main Counting algorithm to look for all steric overlaps *********/
/**************************************************************************/
/**************************************************************************/
double New_Craig_Counting(Mol *M, double eps, unsigned mode)
{
Over *over[MAX_OVER];
char multiple[10][MAX_OVER];
double solid=0.0;
char line[LINELEN];
Atms *at=NULL;
int i=0;
if(M->multi>MAX_OVER) M->multi=MAX_OVER;
for(i=0;imulti;i++)
{
over[i]=NULL;
switch(i+1)
{
case 2 : strcpy(multiple[i],"double"); break;
case 3 : strcpy(multiple[i],"triple"); break;
case 4 : strcpy(multiple[i],"quadruple"); break;
case 5 : strcpy(multiple[i],"quintuple"); break;
case 6 : strcpy(multiple[i],"sextuple"); break;
case 7 : strcpy(multiple[i],"septuple"); break;
case 8 : strcpy(multiple[i],"octuple"); break;
case 9 : strcpy(multiple[i],"nonuple"); break;
case 10 : strcpy(multiple[i],"dectuple"); break;
default : sprintf(multiple[i],"order=%d",i+1); break;
}
}
M->atoms=First_Atom(M->atoms); /* start main loops at first atom */
/****************** ignore all fully overlapped atoms *********************/
for(at=M->atoms;at!=NULL;at=at->next)
{
if((at->stat&MAIN_BIT)&&(Outside_All_Atoms(at)))
at->stat=at->stat|CALC_BIT; else at->stat=at->stat&(FULL_BIT^CALC_BIT);
}
/****************** add all single atom solid angles **********************/
if(mode&VIS_BIT) Out_Message("Calculating all single atom solid angles",O_NEWLN);
solid=Get_All_Single_Solids(M,mode);
/****************** remove all double atom overlaps ***********************/
if(mode&VIS_BIT) Out_Message("Removing all double overlaps",O_NEWLN);
solid=Remove_All_Double_Overlaps(M,&over[1],eps,mode,solid);
/****************** add/subtract all higher order overlaps ****************/
for(i=2;imulti;i++)
{
if(mode&VIS_BIT)
{
if(i/2*2==i) sprintf(line,"Adding all %s overlaps",multiple[i]);
else sprintf(line,"Removing all %s overlaps",multiple[i]);
Out_Message(line,O_NEWLN);
}
solid=Multiple_Overlap(over[i-1],&over[i],over[1]
,i,0,eps,mode,multiple[i],solid);
if(over[i]==NULL) break;
}
/**************************** finish off **********************************/
for(i=1;imulti;i++) Close_All_Overlaps(over[i],multiple[i],mode);
if (solid>4*PI) { Error_Message(E_STOBIG,"New Craig Counting"); solid=4*PI; }
if (solid<0.0) { Error_Message(E_STOSML,"New Craig Counting"); solid=0.0; }
return(solid);
}
/**************************************************************************/
/****************** The End ... *****************************************/
/**************************************************************************/
|