//  hex2.cpp  Aug 19 2002  cleaned up
//  hex2c5   why 66 but no 33 for k6pt???
//      got rid of /r01 factor... - works fine now!
//  hex2c4.cpp  Aug 18  6pt forces cleaned up -- no change
//  hex2c2.cpp  Aug 18  6pt x plane force
//  hex2c1 checking for error (in k3pt) but didn't find any
//  hex2c Aug 17  calc all ksp's for isotropic material
//  hex2b  Aug 17  CdS and CdSe  (don't need same c/a ratio!)
//     Note: uca=1, ucc=0.93
//  hex2  Aug 17
//  hex1f.cpp  Aug 16  display and test multiply matrix inverse.
//    This reveals the problem that ksp1 and ksp2 give exactly
//    the same values of the Cij.  (which makes sense)
//  hex1e  (invert matrix)   just display full Cij matrix
// hex1d.cpp  shorten, same as before
// hex1c.cpp Aug 16  6-point force term
// hex1b.cpp  Aug 16 error fixed; components consistent now
// i.e.  C44=C55  C12=C21  C13=C31 etc.
// hex1.cpp  Aug 16 2002  - use conv. strain def'n
//       (error:  components are not consistent)

void invmat(void);
#define NR 5
double Amat[NR][NR],augmat[NR][NR+NR];
double invAmat[NR][NR];
double testmat[NR][NR];

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

double uca,ucc; // hexag lattice parameters
double blv[3][3]; // 3 Bravais lattice vectors

double dism[3][3];  // space distortion matrix

double e[3][3];  // strain matrix
int ispring,icoef;
double calcc(int ic1,int ic2,int ispring);

double isomat[5][3];

void main(void)
{ // main
clrscr();

uca=1; ucc=0.93;  // Hexagonal lattice constants 'a' and 'c'
//uca=4.30; ucc=7.02;  // CdSe  (Kittel pg. 33) Hexagonal lattice constants 'a' and 'c'
//uca=4.13; ucc=6.75;  // CdS  (Kittel pg. 33) Hexagonal lattice constants 'a' and 'c'

for (ispring=0;ispring<5;ispring++)
for (icoef=0;icoef<5;icoef++)
{ //
Amat[icoef][ispring]=0;
} //

int i1a,i2a;
int straintype;
/*
gotoxy(30,1);cout<<"Options:";
gotoxy(30,2);cout<<"1. xx strain";
gotoxy(30,3);cout<<"2. yy strain";
gotoxy(30,4);cout<<"3. zz strain";
gotoxy(30,5);cout<<"4. yz strain";
gotoxy(30,6);cout<<"5. xz strain";
gotoxy(30,7);cout<<"6. xy strain";
*/
//gotoxy(1,2);cout<<"ispr,ic1,ic2 ? ";
//gotoxy(1,2);cout<<"ispr ? ";
int ispr,ic1,ic2;
//cin>>ispr;
//cin>>ic1;
//if (ic1<1) return;
//if (ic1>6) return;
//cin>>ic2;
cout.precision(5);
double d1;

gotoxy( 8,2);cout<<"xx";
gotoxy(18,2);cout<<"yy";
gotoxy(28,2);cout<<"zz";
gotoxy(38,2);cout<<"yz";
gotoxy(48,2);cout<<"xz";
gotoxy(58,2);cout<<"xy";
gotoxy(1,3);cout<<"xx";
gotoxy(1,4);cout<<"yy";
gotoxy(1,5);cout<<"zz";
gotoxy(1,6);cout<<"yz";
gotoxy(1,7);cout<<"xz";
gotoxy(1,8);cout<<"xy";
gotoxy(1,1);cout<<"Cij for k3pt";
ispr=2;
for (ic1=1;ic1<=6;ic1++)
for (ic2=1;ic2<=6;ic2++)
{ //
d1=calcc(ic1,ic2,ispr);
gotoxy(-4+10*ic2,2+ic1);cout<<d1<<"  ";
} //

gotoxy(1,10);cout<<"Cij for k6pt";
gotoxy( 8,11);cout<<"xx";
gotoxy(18,11);cout<<"yy";
gotoxy(28,11);cout<<"zz";
gotoxy(38,11);cout<<"yz";
gotoxy(48,11);cout<<"xz";
gotoxy(58,11);cout<<"xy";
gotoxy(1,12);cout<<"xx";
gotoxy(1,13);cout<<"yy";
gotoxy(1,14);cout<<"zz";
gotoxy(1,15);cout<<"yz";
gotoxy(1,16);cout<<"xz";
gotoxy(1,17);cout<<"xy";
ispr=0;
for (ic1=1;ic1<=6;ic1++)
for (ic2=1;ic2<=6;ic2++)
{ //
d1=calcc(ic1,ic2,ispr);
gotoxy(-4+10*ic2,11+ic1);cout<<d1<<"  ";
} //

gotoxy(1,18);cout<<"Cij for ksp4";
ispr=4;
for (ic1=1;ic1<=6;ic1++)
for (ic2=1;ic2<=6;ic2++)
{ //
d1=calcc(ic1,ic2,ispr);
gotoxy(-4+10*ic2,19+ic1);cout<<d1<<"  ";
} //

gotoxy(72,20);cout<<"Press";
gotoxy(72,21);cout<<"a key";
gotoxy(72,22);cout<<"to   ";
gotoxy(72,23);cout<<"cont.";
getch();
clrscr();


for (ic1=0;ic1<5;ic1++)
for (ispr=0;ispr<5;ispr++)
{ //.
if (ic1==0) d1=calcc(1,1,ispr);
if (ic1==1) d1=calcc(1,2,ispr);
if (ic1==2) d1=calcc(1,3,ispr);
if (ic1==3) d1=calcc(3,3,ispr);
if (ic1==4) d1=calcc(4,4,ispr);
Amat[ic1][ispr]=d1;
gotoxy(8+10*ispr,2+ic1);cout<<d1<<"  ";
} //.
gotoxy(1,2 );cout<<"C11";
gotoxy(1,3 );cout<<"C12";
gotoxy(1,4 );cout<<"C13";
gotoxy(1,5 );cout<<"C33";
gotoxy(1,6 );cout<<"C44";

invmat();

for (ic1=0;ic1<5;ic1++)
for (ispr=0;ispr<5;ispr++)
{ //.
d1=invAmat[ispr][ic1];
if (fabs(d1)<2e-7) d1=0;
if (fabs(d1)<1e4)
{ //
gotoxy(8+13*ic1,11+ispr);cout<<d1<<"  ";
} //
else
{ //
long int li2=d1;
gotoxy(8+13*ic1,11+ispr);cout<<li2<<"  ";
} //
} //.
gotoxy(10,10);cout<<"C11";
gotoxy(10+13,10);cout<<"C12";
gotoxy(10+2*13,10);cout<<"C13";
gotoxy(10+3*13,10);cout<<"C33";
gotoxy(10+4*13,10);cout<<"C44";
gotoxy(1,11);cout<<"k6pt";
gotoxy(1,12);cout<<"ksp1";
gotoxy(1,13);cout<<"k3pt";
gotoxy(1,14);cout<<"ksp3";
gotoxy(1,15);cout<<"ksp4";

int ic3;

gotoxy(1,18);cout<<"Test: multiply matrix by inverse to get unit matrix:";
for (ic1=0;ic1<=4;ic1++)
for (ic2=0;ic2<=4;ic2++)
{ //..
testmat[ic1][ic2]=0;
for (ic3=0;ic3<=4;ic3++)
{//.
testmat[ic1][ic2]+=Amat[ic1][ic3]*invAmat[ic3][ic2];
}//.
double d1;
d1=testmat[ic1][ic2];
if (fabs(d1)<1e-14) d1=0;
gotoxy(2+15*ic1,19+ic2);cout<<d1<<"  ";
} //..

gotoxy(72,20);cout<<"Press";
gotoxy(72,21);cout<<"a key";
gotoxy(72,22);cout<<"to   ";
gotoxy(72,23);cout<<"cont.";
getch();
clrscr();

gotoxy(1,17);cout<<"Isotropic elastic solid:             limit for                ";
gotoxy(1,18);cout<<"          C11           C44        Poisson ratio              ";
gotoxy(1,19);cout<<"                                                               ";
gotoxy(1,20);cout<<"                                                               ";
gotoxy(1,21);cout<<"                                                               ";
gotoxy(1,22);cout<<"                                                               ";
gotoxy(1,23);cout<<"                                                               ";

for (ic1=0;ic1<2;ic1++)
for (ispr=0;ispr<5;ispr++)
{ //.
isomat[ispr][ic1]=0;
} //.

for (ispr=0;ispr<5;ispr++)
{ //.
isomat[ispr][0]+=   invAmat[ispr][0]; // C11
isomat[ispr][0]+=   invAmat[ispr][1]; // C12
isomat[ispr][1]-= 2*invAmat[ispr][1]; // C12
isomat[ispr][0]+=   invAmat[ispr][2]; // C13
isomat[ispr][1]-= 2*invAmat[ispr][2]; // C13
isomat[ispr][0]+=   invAmat[ispr][3]; // C33
isomat[ispr][1]+=   invAmat[ispr][4]; // C44
isomat[ispr][2]=0.5*(isomat[ispr][1]+2*isomat[ispr][0])/(isomat[ispr][1]+isomat[ispr][0]);
} //.

for (ic1=0;ic1<3;ic1++)
for (ispr=0;ispr<5;ispr++)
{ //.
d1=isomat[ispr][ic1];
if (fabs(d1)<2e-6) d1=0;
gotoxy(9+15*ic1,19+ispr);cout<<d1<<"  ";
} //.
gotoxy(1,19);cout<<"k6pt";
gotoxy(1,20);cout<<"ksp1";
gotoxy(1,21);cout<<"k3pt";
gotoxy(1,22);cout<<"ksp3";
gotoxy(1,23);cout<<"ksp4";

gotoxy(1,11);cout<<"k6pt";
gotoxy(1,12);cout<<"ksp1";
gotoxy(1,13);cout<<"k3pt";
gotoxy(1,14);cout<<"ksp3";
gotoxy(1,15);cout<<"ksp4";
gotoxy(9,10);cout<<"CdSe:";
for (ispr=0;ispr<5;ispr++)
{ //.
double kspi;
kspi=0;
kspi+=invAmat[ispr][0]*74.6e9; // C11 CdSe  E. Rabani J Chem Phys
kspi+=invAmat[ispr][1]*46.1e9; // C12 CdSe  E. Rabani J Chem Phys
kspi+=invAmat[ispr][2]*39.4e9; // C13 CdSe  E. Rabani J Chem Phys
kspi+=invAmat[ispr][3]*81.7e9; // C33 CdSe  E. Rabani J Chem Phys
kspi+=invAmat[ispr][4]*13.0e9; // C44 CdSe  E. Rabani J Chem Phys
gotoxy(8,11+ispr);cout<<kspi<<" ";
} //.

gotoxy(21,10);cout<<"     CdS:";
for (ispr=0;ispr<5;ispr++)
{ //.
double kspi;
kspi=0;
kspi+=invAmat[ispr][0]*90.7e9; // C11 CdS  M. Conry
kspi+=invAmat[ispr][1]*58.1e9; // C12 CdS
kspi+=invAmat[ispr][2]*51.0e9; // C13 CdS
kspi+=invAmat[ispr][3]*93.8e9; // C33 CdS
kspi+=invAmat[ispr][4]*15.04e9; // C44 CdS
gotoxy(25,11+ispr);cout<<kspi<<" ";
} //.


double nuPoisson;
nuPoisson=0.35;
double C11,C44;
C11=90.0e9;
C44=C11*0.5*(1-2*nuPoisson)/(1-nuPoisson);

gotoxy(50,18);cout<<"nu="<<nuPoisson<<" ";
for (ispr=0;ispr<5;ispr++)
{ //.
float ksp=0;
ksp=ksp+isomat[ispr][0]*C11+isomat[ispr][1]*C44;
gotoxy(50,19+ispr);cout<<ksp<<"  ";
} //.

//d1=invAmat[ispr][ic1];
//} //.
getch();

} // main

void invmat(void)
{ // invmat                                                invmat
int i1,i2,i3; double r1;

for (i1=0;i1<NR;i1++)
for (i2=0;i2<NR;i2++)
{ // ..
augmat[i1][i2]=Amat[i1][i2]; augmat[i1][i2+NR]=0;
} // ..
for (i1=0;i1<NR;i1++) augmat[i1][i1+NR]=1;

for (i3=0;i3<NR;i3++)
{ // i3
r1=augmat[i3][i3];
for (i2=0;i2<NR+NR;i2++) augmat[i3][i2]=augmat[i3][i2]/r1;
for (i1=0;i1<NR;i1++)
{ // i1
if (i1!=i3)
{ // not equal
r1=augmat[i1][i3]/augmat[i3][i3];
for (i2=0;i2<NR+NR;i2++) augmat[i1][i2]=augmat[i1][i2]-r1*augmat[i3][i2];
} // not equal
} // i1
} // i3

for (i1=0;i1<NR;i1++)
for (i2=0;i2<NR;i2++)
{ // ...
invAmat[i1][i2]=augmat[i1][NR+i2];
} // ...

return;
} // invmat
//----------------------------------------------------------------------
double calcc(int ic1,int ic2,int ispring2)
{ //   calcc
double cresult;
int straintype=ic2;

double de=0.000001;
int i1a,i2a;
for (i1a=0;i1a<3;i1a++)
for (i2a=0;i2a<3;i2a++)
{ //.
dism[i1a][i2a]=0;
} //.
if (straintype==1) dism[0][0]=de;
if (straintype==2) dism[1][1]=de;
if (straintype==3) dism[2][2]=de;
if (straintype==4) {dism[1][2]=de;dism[2][1]=de;}
if (straintype==5) {dism[0][2]=de;dism[2][0]=de;}
if (straintype==6) {dism[0][1]=de;dism[1][0]=de;}

// Define Bravais lattice vectors
blv[0][0]=uca; blv[0][1]=0; blv[0][2]=0;
blv[1][0]=0.5*uca;
blv[1][1]=0.5*sqrt(3.0)*uca;
blv[2][2]=0;
blv[2][0]=0; blv[2][1]=0; blv[2][2]=ucc;

int ib0,ib1,ib2;
int i0b0,i0b1,i0b2;
double x1,y1,z1;
double x2,y2,z2;
double x01,y01,z01;
double x02,y02,z02;
int bx0,by0; bx0=320; by0=240; double sfx=25;


//gotoxy(1,6);cout<<"Strain tensor (e):";
for (i1a=0;i1a<3;i1a++)
{ //.
e[i1a][i1a]=dism[i1a][i1a];
//gotoxy(1+9*i1a,7+i1a);cout<<e[i1a][i1a]<<" ";
} //.

for (i1a=0;i1a<3;i1a++)
for (i2a=0;i2a<3;i2a++)
{ //.
if (i1a!=i2a)
{ //..
// Most common (inelegant) convention: No factor of a half here:
e[i1a][i2a]=dism[i1a][i2a]+dism[i2a][i1a];
//gotoxy(1+9*i1a,7+i2a);cout<<e[i1a][i2a]<<" ";
} //..
} //.
/*
gotoxy(1,11);cout<<"exx = e1 = "<<e[0][0]<<"  ";
gotoxy(1,12);cout<<"eyy = e2 = "<<e[1][1]<<"  ";
gotoxy(1,13);cout<<"ezz = e3 = "<<e[2][2]<<"  ";
gotoxy(1,14);cout<<"eyz = e4 = "<<e[1][2]<<"  ";
gotoxy(1,15);cout<<"exz = e5 = "<<e[0][2]<<"  ";
gotoxy(1,16);cout<<"exy = e6 = "<<e[0][1]<<"  ";
*/

double mage=1;
int onestraincomponent=0;
if ((onestraincomponent==0)&&(e[0][0]!=0)) {onestraincomponent=1;mage=e[0][0];}
if ((onestraincomponent> 0)&&(e[1][1]!=0)) {onestraincomponent=9;}
if ((onestraincomponent==0)&&(e[1][1]!=0)) {onestraincomponent=2;mage=e[1][1];}
if ((onestraincomponent> 0)&&(e[2][2]!=0)) {onestraincomponent=9;}
if ((onestraincomponent==0)&&(e[2][2]!=0)) {onestraincomponent=3;mage=e[2][2];}
if ((onestraincomponent> 0)&&(e[1][2]!=0)) {onestraincomponent=9;}
if ((onestraincomponent==0)&&(e[1][2]!=0)) {onestraincomponent=4;mage=e[1][2];}
if ((onestraincomponent> 0)&&(e[0][2]!=0)) {onestraincomponent=9;}
if ((onestraincomponent==0)&&(e[0][2]!=0)) {onestraincomponent=5;mage=e[0][2];}
if ((onestraincomponent> 0)&&(e[0][1]!=0)) {onestraincomponent=9;}
if ((onestraincomponent==0)&&(e[0][1]!=0)) {onestraincomponent=6;mage=e[0][1];}

// Find total force on atom 0 0 0
double Fxr,Fyr,Fzr;
double Fxpy,Fypy,Fzpy;
double Fxpz,Fypz,Fzpz;
Fxr=0;
Fyr=0;
Fzr=0;
Fxpy=0;
Fypy=0;
Fzpy=0;
Fxpz=0;
Fypz=0;
Fzpz=0;

double ksp1=0;
double ksp2=0;
double ksp3=0;
double ksp4=0;
double k6pt=0;
double k3pt=0;
double kspi;
if (ispring2==0) k6pt=1;
if (ispring2==1) ksp1=1;
//if (ispring2==2) ksp2=1;
if (ispring2==2) k3pt=1;
if (ispring2==3) ksp3=1;
if (ispring2==4) ksp4=1;

int counter3=0;
double equidist;
for (ispring=1;ispring<=4;ispring++)
{ // ispring=1,2,3,4
if (ispring==1) {equidist=1.0*uca;kspi=ksp1;} // First neighbours
if (ispring==2) {equidist=99*uca;kspi=0;} //
if (ispring==3) {equidist=ucc;kspi=ksp3;} // 3rd neighbours
if (ispring==4) {equidist=sqrt(ucc*ucc+uca*uca);kspi=ksp4;} // 4th neighbours
// Find Force  neighbours
for (i0b0=-2;i0b0<=0;i0b0++)
for (i0b1=-1;i0b1<=0;i0b1++)
for (i0b2=0;i0b2<=0;i0b2++)
for (ib0=-3;ib0<=3;ib0++)
for (ib1=-3;ib1<=3;ib1++)
for (ib2=-1;ib2<=1;ib2++)
{ // ib0 ib1 ib2 -
x1= ib0 *blv[0][0] +ib1*blv[1][0] +ib2*blv[2][0];
y1= ib0 *blv[0][1] +ib1*blv[1][1] +ib2*blv[2][1];
z1= ib0 *blv[0][2] +ib1*blv[1][2] +ib2*blv[2][2];
x01=i0b0*blv[0][0]+i0b1*blv[1][0]+i0b2*blv[2][0];
y01=i0b0*blv[0][1]+i0b1*blv[1][1]+i0b2*blv[2][1];
z01=i0b0*blv[0][2]+i0b1*blv[1][2]+i0b2*blv[2][2];
double r2; r2=(x1-x01)*(x1-x01)+(y1-y01)*(y1-y01)+(z1-z01)*(z1-z01);
int flag2=0;
if ((ispring==1)&&(ib2==i0b2)) flag2=1;
if ((ispring==2)&&(ib2==i0b2)) flag2=1;
if ((ispring==3)&&(i0b0==ib0)&&(i0b1==ib1)) flag2=1;
if ((ispring==4)&&(ib2*ib2==1)&&(i0b2==0)&&(i0b0==0)) flag2=1;
if ((fabs(r2-equidist*equidist)<0.05*equidist*equidist)&&(flag2==1))
{ //..
x2=x1+dism[0][0]*x1+dism[0][1]*y1+dism[0][2]*z1;
y2=y1+dism[1][0]*x1+dism[1][1]*y1+dism[1][2]*z1;
z2=z1+dism[2][0]*x1+dism[2][1]*y1+dism[2][2]*z1;
x02=x01+dism[0][0]*x01+dism[0][1]*y01+dism[0][2]*z01;
y02=y01+dism[1][0]*x01+dism[1][1]*y01+dism[1][2]*z01;
z02=z01+dism[2][0]*x01+dism[2][1]*y01+dism[2][2]*z01;
double r02; r02=(x2-x02)*(x2-x02)+(y2-y02)*(y2-y02)+(z2-z02)*(z2-z02);
double r01;r01=sqrt(r02);
double Fmag=kspi*(r01-equidist);
double F01x,F01y,F01z;
F01x=Fmag*(x2-x02)/r01;
F01y=Fmag*(y2-y02)/r01;
F01z=Fmag*(z2-z02)/r01;
if ((x1>0)&&(x01<=0)&&(fabs(y01)<0.9*uca)&&(y01<=0))
{ //.
//if ((ispring==4))
//{ //..
//counter3++;
//gotoxy(50,counter3);cout<<"xxx "<<ib0<<" "<<ib1<<" "<<ib2<<" "<<i0b0<<" "<<i0b1<<" "<<i0b2<<" ";
//} //..
  Fxr+=F01x;
  Fyr+=F01y;
  Fzr+=F01z;
} //.
if ((y1>0)&&(y01<=0)&&(fabs(x01)<0.9*uca)&&(x01<=0))
{ //.
  Fxpy+=F01x;
  Fypy+=F01y;
  Fzpy+=F01z;
} //.
if ((z1>0)&&(z01<=0)&&(fabs(x01)<0.4*uca)&&(fabs(y01)<0.4*uca))
{ //.
  Fxpz+=F01x;
  Fypz+=F01y;
  Fzpz+=F01z;
} //.
} //..
} // ib0 ib1 ib2
} // ispring=1,2,3,4

int iway;

// Enumerate all 6-atom clusters:
for (ib0=-3;ib0<=3;ib0++)
for (ib1=-3;ib1<=3;ib1++)
for (ib2=-3;ib2<=3;ib2++)
for (iway=-1;iway<=1;iway=iway+2)
{ // ib0 ib1 ib2 - 6 point force - k6pt
// First, locate center of mass of the cluster:
double cmx,cmy,cmz,xmin,xmax,ymin,ymax,zmin,zmax;
cmx=0;
cmy=0;
cmz=0;
xmin=100;xmax=-100;
ymin=100;ymax=-100;
zmin=100;zmax=-100;
int ib0a,ib1a,ib2a,ipt;
for (ib2a=ib2;ib2a<=ib2+1;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
x2=x1+dism[0][0]*x1+dism[0][1]*y1+dism[0][2]*z1;
y2=y1+dism[1][0]*x1+dism[1][1]*y1+dism[1][2]*z1;
z2=z1+dism[2][0]*x1+dism[2][1]*y1+dism[2][2]*z1;
cmx=cmx+x2;cmy=cmy+y2;cmz=cmz+z2;
if (x1>xmax) xmax=x1;
if (y1>ymax) ymax=y1;
if (z1>zmax) zmax=z1;
if (x1<xmin) xmin=x1;
if (y1<ymin) ymin=y1;
if (z1<zmin) zmin=z1;
} // ib2a
cmx=cmx/6.0;
cmy=cmy/6.0;
cmz=cmz/6.0;
//putpixel(bx0+sfx*cmx,by0-sfx*cmy,10+iway);

double sumr2=0,r2;
for (ib2a=ib2;ib2a<=ib2+1;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
x2=x1+dism[0][0]*x1+dism[0][1]*y1+dism[0][2]*z1;
y2=y1+dism[1][0]*x1+dism[1][1]*y1+dism[1][2]*z1;
z2=z1+dism[2][0]*x1+dism[2][1]*y1+dism[2][2]*z1;
r2=(cmx-x2)*(cmx-x2)+(cmy-y2)*(cmy-y2)+(cmz-z2)*(cmz-z2);
sumr2=sumr2+r2;
} // ib2a, ipt

double Fmag;
double F01x,F01y,F01z,r01;

sumr2=sumr2-6.0*((uca*uca/3.0)+(ucc*ucc/4.0));
Fmag=sumr2*k6pt;

for (ib2a=ib2;ib2a<=ib2+1;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
x2=x1+dism[0][0]*x1+dism[0][1]*y1+dism[0][2]*z1;
y2=y1+dism[1][0]*x1+dism[1][1]*y1+dism[1][2]*z1;
z2=z1+dism[2][0]*x1+dism[2][1]*y1+dism[2][2]*z1;
r01=sqrt((x2-cmx)*(x2-cmx)+(y2-cmy)*(y2-cmy)+(z2-cmz)*(z2-cmz));
//F01x=Fmag*(cmx-x2)/r01;
//F01y=Fmag*(cmy-y2)/r01;
//F01z=Fmag*(cmz-z2)/r01;
F01x=Fmag*(cmx-x2); // fixed Aug 18
F01y=Fmag*(cmy-y2); // fixed Aug 18
F01z=Fmag*(cmz-z2); // fixed Aug 18

if ((xmin<=0)&&(xmax>0)&&(ib2a==0)&&(y1<=0)&&(y1>=-0.9*uca)&&(x1<=0))
{ //.. changed in hex2c2.cpp
Fxr+=F01x;
Fyr+=F01y;
Fzr+=F01z;
} //..

if ((ymin<=0)&&(ymax>0)&&(ib2a==0)&&(ib0a==0)&&(y1<=0))
{ //   changed in hex2c3.cpp
Fxpy+=F01x;
Fypy+=F01y;
Fzpy+=F01z;
} //...

if ((zmin<=0)&&(zmax>0)&&(ib2a==0)&&(ib1a==0)&&(ib0a==0))
{ //.. changed in hex2c4.cpp
Fxpz+=F01x;
Fypz+=F01y;
Fzpz+=F01z;
} //..

} // ib2a, ipt


} // ib0 ib1 ib2 - 6 point force - k6pt

int count333=0;
// Enumerate all 3-atom clusters:
for (ib0=-3;ib0<=3;ib0++)
for (ib1=-3;ib1<=3;ib1++)
for (ib2=-3;ib2<=3;ib2++)
for (iway=-1;iway<=1;iway=iway+2)
{ // ib0 ib1 ib2 - 3 point force - k3pt
// First, locate center of mass of the cluster:
// Aug 18: also locate undistorted max and min of cluster
double cmx,cmy,cmz,xmin,xmax,ymin,ymax,zmin,zmax;
cmx=0;
cmy=0;
cmz=0;
xmin=100;xmax=-100;
ymin=100;ymax=-100;
zmin=100;zmax=-100;
int ib0a,ib1a,ib2a,ipt;
for (ib2a=ib2;ib2a<=ib2+0;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
x2=x1+dism[0][0]*x1+dism[0][1]*y1+dism[0][2]*z1;
y2=y1+dism[1][0]*x1+dism[1][1]*y1+dism[1][2]*z1;
z2=z1+dism[2][0]*x1+dism[2][1]*y1+dism[2][2]*z1;
cmx=cmx+x2;cmy=cmy+y2;cmz=cmz+z2;
if (x1>xmax) xmax=x1;
if (y1>ymax) ymax=y1;
if (z1>zmax) zmax=z1;
if (x1<xmin) xmin=x1;
if (y1<ymin) ymin=y1;
if (z1<zmin) zmin=z1;
} // ib2a
cmx=cmx/3.0;
cmy=cmy/3.0;
cmz=cmz/3.0;
//putpixel(bx0+sfx*cmx,by0-sfx*cmy,10+iway);

double sumr2=0,r2;
for (ib2a=ib2;ib2a<=ib2+0;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
x2=x1+dism[0][0]*x1+dism[0][1]*y1+dism[0][2]*z1;
y2=y1+dism[1][0]*x1+dism[1][1]*y1+dism[1][2]*z1;
z2=z1+dism[2][0]*x1+dism[2][1]*y1+dism[2][2]*z1;
r2=(cmx-x2)*(cmx-x2)+(cmy-y2)*(cmy-y2)+(cmz-z2)*(cmz-z2);
sumr2=sumr2+r2;
} // ib2a, ipt

double Fmag;
double F01x,F01y,F01z,r01;

sumr2=sumr2-3.0*(uca*uca/3.0);
Fmag=sumr2*k3pt;

for (ib2a=ib2;ib2a<=ib2+0;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
x2=x1+dism[0][0]*x1+dism[0][1]*y1+dism[0][2]*z1;
y2=y1+dism[1][0]*x1+dism[1][1]*y1+dism[1][2]*z1;
z2=z1+dism[2][0]*x1+dism[2][1]*y1+dism[2][2]*z1;
r01=sqrt((x2-cmx)*(x2-cmx)+(y2-cmy)*(y2-cmy)+(z2-cmz)*(z2-cmz));
//F01x=Fmag*(cmx-x2)/r01;
//F01y=Fmag*(cmy-y2)/r01;
//F01z=Fmag*(cmz-z2)/r01;
F01x=Fmag*(cmx-x2); // fixed Aug 18
F01y=Fmag*(cmy-y2); // fixed Aug 18
F01z=Fmag*(cmz-z2); // fixed Aug 18
// Aug 18: It should be force on 000 and other inside atoms,
// depending which surface we are looking at:

if ((ib2a==0)&&(xmax>0)&&(xmin<=0)&&(x1<=0)&&(y1<=0)&&(y1>=-0.9*uca))
{ //.
  Fxr+=F01x;
  Fyr+=F01y;
  Fzr+=F01z;
} //.

if ((ib0a==0)&&(ib1a==0)&&(ib2a==0))
{ // force on 000
if (cmy>0.0)
{ //...
  Fxpy+=F01x;
  Fypy+=F01y;
  Fzpy+=F01z;
} //...
if (cmz>0.0)
{ //..
  Fxpz+=F01x; // There are no z-forces anyways, I think??
  Fypz+=F01y;
  Fzpz+=F01z;
} //..
} // force on 000

} // ib2a, ipt


} // ib0 ib1 ib2 - 3 point force - k3pt


double sigxx,sigyy,sigzz;
double sigxy,sigyz,sigxz;
double sigyx,sigzy,sigzx;

sigxx=Fxr/(ucc*sqrt(3.0)*uca);
sigxy=Fyr/(ucc*sqrt(3.0)*uca);
sigxz=Fzr/(ucc*sqrt(3.0)*uca);
sigyx=Fxpy/(ucc*uca);
sigyy=Fypy/(ucc*uca);
sigyz=Fzpy/(ucc*uca);
sigzx=Fxpz/(uca*uca*sqrt(3.0)*0.5);
sigzy=Fypz/(uca*uca*sqrt(3.0)*0.5);
sigzz=Fzpz/(uca*uca*sqrt(3.0)*0.5);

double C1n=sigxx/mage; if (fabs(C1n)<1e-5) C1n=0;
double C2n=sigyy/mage; if (fabs(C2n)<1e-5) C2n=0;
double C3n=sigzz/mage; if (fabs(C3n)<1e-5) C3n=0;
double C4n=sigyz/mage; if (fabs(C4n)<1e-5) C4n=0;
double C5n=sigxz/mage; if (fabs(C5n)<1e-5) C5n=0;
double C6n=sigxy/mage; if (fabs(C6n)<1e-5) C6n=0;
if (ic1==1) cresult=C1n;
if (ic1==2) cresult=C2n;
if (ic1==3) cresult=C3n;
if (ic1==4) cresult=C4n;
if (ic1==5) cresult=C5n;
if (ic1==6) cresult=C6n;
/*
gotoxy(20,18);cout<<"C1"<<onestraincomponent<<" = "<<C1n<<" ";
gotoxy(20,19);cout<<"C2"<<onestraincomponent<<" = "<<C2n<<" ";
gotoxy(20,20);cout<<"C3"<<onestraincomponent<<" = "<<C3n<<" ";
gotoxy(20,21);cout<<"C4"<<onestraincomponent<<" = "<<C4n<<" ";
gotoxy(20,22);cout<<"C5"<<onestraincomponent<<" = "<<C5n<<" ";
gotoxy(20,23);cout<<"C6"<<onestraincomponent<<" = "<<C6n<<" ";
*/

return cresult;
} // calcc