// same8ptf.cpp  Aug 5, 2002
//   check that 8 point force comes from 8 point potential

#include <iostream.h>
#include <conio.h>
#include <math.h>
#include <stdlib.h>
//#include <.h>
//#include <.h>

double k8pt=1.5;
double uca=1.3;
double r[8][3];

void U(void);
void F(int ij3,int ip3);
double Uret,Fret;
double RCM[3];
double D;

void main(void)
{ //
clrscr();
randomize();
int ij,ip;
for (ij=0;ij<8;ij++)
{ // ij
for (ip=0;ip<3;ip++)
{ // ip
r[ij][ip]=0.0001*uca*random(10000);
} // ip
} // ij

double U000;
U();
U000=Uret;
double dx=0.01;
double diff5;
double Unew,Uplus,Uminus;   long int lidiff;
int ij3,ip3;

LB1:
gotoxy(40,1);cout<<"dx = "<<dx<<"  ";

gotoxy(1,3);cout<<"Force from force formula:";
for (ij3=0;ij3<8;ij3++)
{ // ij3
gotoxy(1,5+ij3);cout<<"F["<<ij3<<"] = ";
for (ip3=0;ip3<3;ip3++)
{ // ip3
F(ij3,ip3);
gotoxy(8+10*ip3,5+ij3);cout<<Fret<<"  ";
} // ip3
} // ij3

gotoxy(40,3);cout<<"Force from potential derivative:";
for (ij3=0;ij3<8;ij3++)
{ // ij3
gotoxy(40,5+ij3);cout<<"F["<<ij3<<"] = ";
for (ip3=0;ip3<3;ip3++)
{ // ip3
r[ij3][ip3]=r[ij3][ip3]+dx;
U(); Uplus=Uret;
r[ij3][ip3]=r[ij3][ip3]-2*dx;
U(); Uminus=Uret;
r[ij3][ip3]=r[ij3][ip3]+dx;
gotoxy(40+7+10*ip3,5+ij3);cout<<-(Uplus-Uminus)/(2*dx)<<"  ";
} // ip3
} // ij3

gotoxy(1,15);cout<<"Difference times 100,000,000";
for (ij3=0;ij3<8;ij3++)
{ // ij3
gotoxy(1,16+ij3);cout<<"d["<<ij3<<"] = ";
for (ip3=0;ip3<3;ip3++)
{ // ip3
r[ij3][ip3]=r[ij3][ip3]+dx;
U(); Uplus=Uret;
r[ij3][ip3]=r[ij3][ip3]-2*dx;
U(); Uminus=Uret;
r[ij3][ip3]=r[ij3][ip3]+dx;
F(ij3,ip3);
diff5=Fret+(Uplus-Uminus)/(2*dx);
diff5=diff5*1e8;
lidiff=diff5;
gotoxy(1+7+15*ip3,16+ij3);cout<<lidiff<<"  ";
} // ip3
} // ij3



gotoxy(50,22);cout<<"Press 'q' to quit";
gotoxy(50,24);cout<<"Press spacebar to repeat";
gotoxy(50,25);cout<<"with dx 20% smaller";
char ch1;
ch1=getch();
if (ch1=='q') return;
dx=dx*0.8;
goto LB1;
} //

void U(void)
{ //  U
double Uans;
int ij,ip;
RCM[0]=0;
RCM[1]=0;
RCM[2]=0;
for (ij=0;ij<8;ij++)
{ // ij
for (ip=0;ip<3;ip++)
{ // ip
RCM[ip]+=r[ij][ip];
} // ip
} // ij
RCM[0]=RCM[0]/8.0;
RCM[1]=RCM[1]/8.0;
RCM[2]=RCM[2]/8.0;
D=0;
for (ij=0;ij<8;ij++)
{ // ij
for (ip=0;ip<3;ip++)
{ // ip
D = D + (RCM[ip]-r[ij][ip])*(RCM[ip]-r[ij][ip]);
} // ip
} // ij
Uans = k8pt/(48.0*uca*uca);
Uans=Uans*(D-6*uca*uca)*(D-6*uca*uca);
Uret=Uans;
} //  U

void F(int ij3,int ip3)
{ //  F
double Fijip;
int ij,ip;
RCM[0]=0;
RCM[1]=0;
RCM[2]=0;
for (ij=0;ij<8;ij++)
{ // ij
for (ip=0;ip<3;ip++)
{ // ip
RCM[ip]+=r[ij][ip];
} // ip
} // ij
RCM[0]=RCM[0]/8.0;
RCM[1]=RCM[1]/8.0;
RCM[2]=RCM[2]/8.0;
D=0;
for (ij=0;ij<8;ij++)
{ // ij
for (ip=0;ip<3;ip++)
{ // ip
D = D + (RCM[ip]-r[ij][ip])*(RCM[ip]-r[ij][ip]);
} // ip
} // ij
Fijip=k8pt/(12.0*uca*uca);
Fijip=Fijip*(D-6.0*uca*uca);
Fijip=Fijip*(RCM[ip3]-r[ij3][ip3]);
Fret=Fijip;
//return Fijip;
} //  F