icorot.f
c icosahedron rotator program 0.1a implicit real*8(a-h,o-z) parameter(n=60) double precision m,m2,xtemp,ytemp,ztemp,theta,a,b,p double precision conv,switch,rad dimension x(n),y(n),z(n) open(12,file=’c60.xyz’) open(13,file=’c60copy.xyz’) open(14,file=’c60ico.xyz’) switch=0.0d0 p=20.90515744515d0 10 format(1x,f24.20,1x,f24.20,1x,f24.20) 81 format(f24.20,1x,f24.20,1x,I10) do i=1,n read(12,10)x(i),y(i),z(i) enddo rad=3.1415926540d0/180.0d0 do i=1,n x(i)=x(i)*rad y(i)=y(i)*rad z(i)=z(i)*rad enddo do i=1,n write(13,10)x(i),y(i),z(i) enddo do i=1,9999999 if (switch.eq.1) goto 91 p=p*rad+1.0E-14*rad a=(y(1)/rad)*sin(p)+(z(1)/rad)*cos(p) b=(y(2)/rad)*sin(p)+(z(2)/rad)*cos(p) conv=abs(a-b) p=p/rad c write(6,*)p,a,b,conv,i write(6,81)p,conv,i if (conv.lt.1.0E-15) switch=1 enddo 91 continue p=p*rad do i=1,n ytemp=0.0d0 ztemp=0.0d0 ytemp=y(i)*cos(p)-z(i)*sin(p) ztemp=y(i)*sin(p)+z(i)*cos(p) y(i)=ytemp z(i)=ztemp enddo write(13,*)’rotated’ do i=1,n x(i)=x(i)/rad y(i)=y(i)/rad z(i)=z(i)/rad enddo do i=1,n write(13,10)x(i),y(i),z(i) write(14,10)x(i),y(i),z(i) enddo stop end |
pent.f
c implicit real*8 (a-h,o-z)
character*4 itex,itexh,itexb,ibas,ibash,istg
integer p, s
double precision aconst,r,tconst,zdf,test,pol,rr,dota,dotb
double precision ratph, rstart, rstep, rpent, rhex,switch
double precision conv
parameter (n=5)
parameter (s=200)
dimension x(12),y(12),z(12),xc(20),yc(20),zc(20)
dimension xic(30),yic(30),zic(30),xp(12),yp(12),op(12)
dimension xb(12),yb(12),zb(12),xbh(12),ybh(12),zbh(12)
dimension xh(20),yh(20),zh(20),xi(60),yi(60),zi(60)
Dimension F(3), alp(3,3), xyz(n,3), b(3), d(3), t(3,3)
dimension af(3), at(3,3), stmew(n,3), cor(n,n,3), old(3)
dimension c(3), e(3), g(3), stold(n,3), pent(5)
dimension vx(n), vy(n), vz(n),stdif(n,3)
data xp/0.5d0,-0.5d0,0.0d0,5.0d0,5.0d0,0.0d0,-5.0d0,0.0d0,
10.5d0,0.0d0,-5.0d0,-0.5d0/
data yp/0.0d0,0.0d0,5.0d0,0.5d0,-0.5d0,-5.0d0,0.5d0,5.0d0,
10.0d0,-5.0d0,-0.5d0,0.0d0/
data op/5.0d0,5.0d0,0.5d0,0.0d0,0.0d0,0.5d0,0.0d0,-0.5d0,
1-5.0d0,-0.5d0,0.0d0,-5.0d0/
data itex/'C '/
data itexh/'H '/
data itexb/'B '/
data ibas/'P '/
data ibash/'S '/
data istg/' 12C'/
data conv/0.5291771d0/
zdf=0.1d0
aconst=0.25d0
m=0.0d0
switch=0.0d0
conv=0.0d0
open(14,file='cplay.xyz')
open(15,file='q.in')
open(13,file='c60.xyz')
open(7,file='f.dat')
open(8,file='data.out')
open(10,file='qrun.out')
open(12,file='debug.dat')
open(16,file='c60.dat')
open(17,file='spread.out')
read(15,*)ratph,rstart,rstep
write(6,*)'enter pent numbers'
read(5,*)pent
c The q-loop!
do q=1,1
do i=1,12
x(i)=xp(i)
y(i)=yp(i)
z(i)=op(i)
enddo
c chris here
rewind 16
rewind 7
if (q.ne.1) rstart=rstart+rstep
rpent=rstart
rhex=ratph*rpent
write(10,*)'ratph',ratph,'rstart',rstart
write(10,*)'rpent',rpent,'rhex',rhex,'rstep',rstep
read(16,*)rbb,rbh
read(16,*)rcc,rch
write(10,*)'rbb',rbb,'rbh',rbh
write(10,*)'rcc',rcc,'rch',rch
c chris ends here
rbb=rbb/conv
rbh=rbh/conv
rcc=rcc/conv
rch=rch/conv
rpent=rpent/conv
rhex=rhex/conv
char=6.0d0
charh=1.0d0
charb=5.0d0
phi=0.5d0*(1.0d0+dsqrt(5.0d0))
i=0
j=0
k=0
do 1 i=1,12
if(x(i).gt.2.0d0)x(i)=phi/2.0d0
if(x(i).lt.-2.0d0)x(i)=-phi/2.0d0
if(y(i).gt.2.0d0)y(i)=phi/2.0d0
if(y(i).lt.-2.0d0)y(i)=-phi/2.0d0
if(z(i).gt.2.0d0)z(i)=phi/2.0d0
if(z(i).lt.-2.0d0)z(i)=-phi/2.0d0
1 continue
c read the faces of the icosahedron
do 2 ic=1,20
read(16,*)i, j, k
xc(ic)=(x(i)+x(j)+x(k))/phi
yc(ic)=(y(i)+y(j)+y(k))/phi
zc(ic)=(z(i)+z(j)+z(k))/phi
2 continue
c write(10,*)'x(1) at 2',x(1)
write(6,89)
89 format(1X,' WARNING!! INPUT ANGSTROM OUTPUT BOHR!!!!! ')
write(6,90)
90 format(1X,' ICOSAHEDRON OF UNIT SIDE ')
do 3 i=1,9
3 write(6,100)itex,i,ibas,x(i),y(i),z(i),char
do 33 i=10,12
33 write(6,101)itex,i,ibas,x(i),y(i),z(i),char
100 format(A1,I1,2X,2X,A2,3X,4F12.9)
101 format(A1,I2,1X,2X,A2,3X,4F12.9)
write(6,80)rbb
80 format(1X,' B12 POSITIONS FOR RBB= ',1X,F10.6)
c scale for boron cage of B12H12
c ze problem is ze ere
do 81 i=1,12
xb(i)=x(i)*rbb
yb(i)=y(i)*rbb
zb(i)=z(i)*rbb
if(i.le.9)then
write(6,100)itexb,i,ibas,xb(i),yb(i),zb(i),charb
else
write(6,101)itexb,i,ibas,xb(i),yb(i),zb(i),charb
endif
81 continue
c scale for H cage
rad=dsqrt(xb(1)**2.0d0+yb(1)**2.0d0+zb(1)**2.0d0)
c JIZZ!! taking off 'd' what will it do??
write(6,82)rbb,rbh
82 format(1X,' H12 POSITIONS FOR RBB,RBH= ',1X,F10.6,1X,F10.6)
do 83 i=1,12
xbh(i)=x(i)*rbb*(rbh+rad)/rad
ybh(i)=y(i)*rbb*(rbh+rad)/rad
zbh(i)=z(i)*rbb*(rbh+rad)/rad
if(i.le.9)then
write(6,100)itexh,i,ibash,xbh(i),ybh(i),zbh(i),charh
else
write(6,101)itexh,i,ibash,xbh(i),ybh(i),zbh(i),charh
endif
83 continue
write(6,91)
91 format(1X,' DODECAHEDRON OF UNIT SIDE ')
do 4 i=1,9
4 write(6,100)itex,i,ibas,xc(i),yc(i),zc(i),char
do 34 i=10,20
34 write(6,101)itex,i,ibas,xc(i),yc(i),zc(i),char
write(6,92)rcc
92 format(1X,' C20 WITH SIDE RCC=',1X,F10.6)
c scale dodecahedron
do 5 i=1,20
xc(i)=rcc*xc(i)
yc(i)=rcc*yc(i)
zc(i)=rcc*zc(i)
if(i.le.9)then
write(6,100)itex,i,ibas,xc(i),yc(i),zc(i),char
else
write(6,101)itex,i,ibas,xc(i),yc(i),zc(i),char
endif
5 continue
rad=dsqrt(xc(1)**2.0d0+yc(1)**2.0d0+zc(1)**2.0d0)
write(6,93)rcc,rch
93 format(1X,' H20 POSITIONS FOR RCC,RCH =',1X,F10.6,1X,F10.6)
do 6 i=1,20
xh(i)=xc(i)*(rch+rad)/rad
yh(i)=yc(i)*(rch+rad)/rad
zh(i)=zc(i)*(rch+rad)/rad
if(i.le.9)then
write(6,100)itexh,i,ibash,xh(i),yh(i),zh(i),charh
else
write(6,101)itexh,i,ibash,xh(i),yh(i),zh(i),charh
endif
6 continue
write(12,94)rpent,rhex
94 format(1X,' C60 POSITIONS FOR RPENT,RHEX = ',2(1X,F10.6))
alpha=rpent/(2.0D0*rpent+rhex)
rbig=rhex+2.0d0*rpent
c scale master icosahedron
do 7 i=1,12
x(i)=x(i)*rbig
y(i)=y(i)*rbig
z(i)=z(i)*rbig
7 continue
c read edges of icosahedron
do 8 i=1,30
read(16,*)j,k
xi(i)=x(j)+alpha*(x(k)-x(j))
yi(i)=y(j)+alpha*(y(k)-y(j))
zi(i)=z(j)+alpha*(z(k)-z(j))
xi(i+30)=x(k)+alpha*(x(j)-x(k))
yi(i+30)=y(k)+alpha*(y(j)-y(k))
zi(i+30)=z(k)+alpha*(z(j)-z(k))
xic(i)=x(j)+x(k)
yic(i)=y(j)+y(k)
zic(i)=z(j)+z(k)
ric=xic(i)*xic(i)+yic(i)*yic(i)+zic(i)*zic(i)
xic(i)=xic(i)/ric
yic(i)=yic(i)/ric
zic(i)=zic(i)/ric
8 continue
c altered...
123 format(4x,d20.16,1x,d20.16,1x,d20.16)
do 888 i=1,9
c 888 write(12,100)itex,i,ibas,xi(i),yi(i),zi(i),char
888 write(14,123)xi(i),yi(i),zi(i)
do 889 i=10,60
c 889 write(12,101)itex,i,ibas,xi(i),yi(i),zi(i),char
889 write(14,123)xi(i),yi(i),zi(i)
c altered...
write(10,*)pent
vx(1)=xi(pent(1))
vy(1)=yi(pent(1))
vz(1)=zi(pent(1))
vx(2)=xi(pent(2))
vy(2)=yi(pent(2))
vz(2)=zi(pent(2))
vx(3)=xi(pent(3))
vy(3)=yi(pent(3))
vz(3)=zi(pent(3))
vx(4)=xi(pent(4))
vy(4)=yi(pent(4))
vz(4)=zi(pent(4))
vx(5)=xi(pent(5))
vy(5)=yi(pent(5))
vz(5)=zi(pent(5))
zero=0.0D0
izero=0
do 977 i=1,60
xi(i)=xi(i)*conv
yi(i)=yi(i)*conv
977 zi(i)=zi(i)*conv
do 988 i=1,60
988 write(6,111)istg,izero,izero,xi(i),yi(i),zi(i),izero,izero,zero
111 format(1X,A4,2I5,3F10.7,2I5,F10.4)
write(6,95)
95 format(1X,' ICOSIDODECAHEDRON OF UNIT RADIUS ')
do 444 i=1,30
444 write(6,445)i,xic(i),yic(i),zic(i)
445 format(1x,i3,3f16.10)
c polar coordinates of c60 in this orientation
write(6,2222)
2222 format(1X,' POLAR COORDINATES OF C60 ')
pi=3.1415926536d0
do 989 i=1,60
r=sqrt(xi(i)*xi(i)+yi(i)*yi(i)+zi(i)*zi(i))
xi(i)=xi(i)/r
yi(i)=yi(i)/r
zi(i)=zi(i)/r
temp=zi(i)
if(temp.gt.1.0d0)temp=1.0d0
if(temp.lt.-1.0d0)temp=-1.0d0
tang=acos(temp)
s1=sin(tang)
test=s1*s1
if(test.lt.1.0d-6)goto 2727
temp=xi(i)/s1
if(temp.gt.1.0d0)temp=1.0d0
if(temp.lt.-1.0d0)temp=-1.0d0
ping=acos(temp)
goto 2828
2727 ping=0.0d0
2828 continue
if(yi(i).lt.0.0d0)ping=2.0d0*pi-ping
tang=tang*180.0d0/pi
ping=ping*180.0d0/pi
write(6,1111)i,xi(i),yi(i),zi(i),tang,ping
1111 format(i3,5f14.4)
989 continue
c chris's program starts here
c *****************************************************
do i=1,3
read(7,*)f(i)
write(6,*)'reading f from f.dat',f(i)
enddo
10 format(4x,d20.16,1x,d20.16,1x,d20.16)
rr=SQRT((vx(1)**2)+(vy(1)**2)+(vz(1)**2))
do i=1,n
xyz(i,1)=vx(i)
xyz(i,2)=vy(i)
xyz(i,3)=vz(i)
enddo
c *****************************************************
c set up alpha
do i=1,3
do j=1,3
if (i.eq.j) then
alp(i,j)=aconst
else
alp(i,j)=0.0d0
endif
enddo
enddo
c ******************************************************
zdff=1.0d0
c h loop numbet of total cycles set by 's'
do h=1,10
if (h.ne.1) zdff=zdf
if(switch.eq.1) goto 67
do l=1,n
do m=1,n
if (m.eq.l) then
goto 98
else
endif
c ******************************************************
c d(3) displacement vector between l and m
do i=1,3
d(i) = 0.0d0
d(i) = xyz(l,i)-xyz(m,i)
enddo
c r distance between i and j
r=SQRT ((d(1)**2)+(d(2)**2)+(d(3)**2))
c af multiplying alp(3,3) by F(3) to get af(3)
do i=1,3
af(i)=0.0d0
do j=1,3
af(i)=af(i)+(alp(i,j)*(f(j)))
enddo
enddo
c t construction
tconst=3.0d0/(r**5)
do i=1,3
do j=1,3
if (i.eq.j) then
t(i,j)=tconst*((D(i)**2)-((1.0d0/3.0d0)*(r**2)))
else
t(i,j)=tconst*(d(i)*(d(j)))
end if
enddo
enddo
c at construction
do i=1,3
do j=1,3
at(i,j)=0.0d0
do k=1,3
at(i,j)=at(i,j)+alp(i,k)*t(k,j)
enddo
enddo
enddo
c *****************************************************
do i=1,3
if (h.eq.1) then
stmew(l,i)=af(i)
old(i)=stmew(l,i)
else
old(i)=stmew(l,i)
endif
enddo
do i=1,3
cor(l,m,i)=0.0d0
do j=1,3
cor(l,m,i)=cor(l,m,i)+(at(i,j)*(old(j)))
enddo
enddo
98 continue
enddo
c write(10,*)'r',r
enddo
c ************* end l and m loops *********************
do l=1,n
do i=1,3
stmew(l,i)=af(i)
if (h.eq.1) stold(l,i)=af(i)
do j=1,n
if (j.eq.l) then
write(10,*)h,conv
else
stmew(l,i)=stmew(l,i)+cor(j,l,i)
write(6,*)stmew(l,i)
endif
enddo
enddo
enddo
do l=1,n
do i=1,3
stmew(l,i)=zdff*stmew(l,i)+(1.0d0-zdff)*stold(l,i)
stdif(l,i)=abs(stmew(l,i)-stold(l,i))
dota=dota+stdif(l,i)*stdif(l,i)
dotb=dotb+stmew(l,i)*stmew(l,i)
conv=dota/dotb
if (conv.lt.0.001d0) switch=1
stold(l,i)=stmew(l,i)
enddo
enddo
enddo
67 continue
c ************ end of h loop **************************
110 do l=1,n
c write(12,*)'particle',l
c do i=1,3
c write(12,*)stmew(l,i)
c enddo
enddo
c finding polarisability
pol=0.0d0
do l=1,n
pol=pol+stmew(l,3)
enddo
write(12,*)'r-displacement from origin',rr
write(12,*)'sum dipoles
|
vsh.f
c VSH 0.9d implicit real*8 (a-h,o-z) parameter (n=60) double precision soconst, v4k, v5k, r, kk, dotp,v6k,v7k,v8k double precision v9k,v100k,v110k,k0,k1,k5,k7,k9,k11,k13 double precision v130k,v3k,k3,v150k,k15,k17 double precision a1,a2,a3,a4,a5,a6,a7,a8,a9 dimension xyz(n,3),x(n),y(n),z(n),so(n,3),v10(n,3),v20(n,3) dimension v30(n,3), v40(n,3), v50(n,3), vnorm(n,3), c60(n,3) dimension vorth(n,3),xi(n),yi(n),zi(n),v60(n,3),v70(n,3),v80(n,3) dimension v90(n,3),v100(n,3),v110(n,3),cage(n,3),v130(n,3) dimension sum(3),v150(n,3),v170(n,3) open(14,file='test.out') c open(13,file='c60ico.xyz') open(13,file='c60.xyz') open(12,file='vsh.out') open(11,file='vsh.in') open(10,file='check.out') open(9,file='vshc60.out') c ********************************************************** c reading in co-ordinates of c60 from c60.xyz do i=1,n read(13,10)x(i), y(i), z(i) enddo c reading in dipoles from vsh.in a1=1.0d0 write(14,28)a1,a1,a1 do i=1,n read(11,28)xi(i),yi(i),zi(i) c60(i,1)=xi(i) c60(i,2)=yi(i) c60(i,3)=zi(i) enddo 10 format(1x,d24.20,1x,d24.20,1x,d24.20) 40 format(2x,f23.20,1x,f23.20,1x,f23.20) 50 format(2x,f23.20,2x,f23.20,2x,f23.20) 25 format(4x,d20.12,1x,d20.12,1x,d20.12) 28 format(1x,d24.16,1x,d24.16,1x,d24.16) r=SQRT ((x(1)**2)+(y(1)**2)+(z(1)**2)) write(6,*)'r',r r=SQRT ((x(2)**2)+(y(2)**2)+(z(2)**2)) write(6,*)'r',r r=SQRT ((x(3)**3)+(y(3)**2)+(z(3)**2)) write(6,*)'r',r do i=1,n x(i)=x(i)/r y(i)=y(i)/r z(i)=z(i)/r enddo do i=1,10 write(6,*)'z',z(i) enddo c ********************************************************** write(12,*)'c60 dipole field' call output(c60,n) c So (0,0,1) write(12,*)'So' c soconst=1.0d0/(60.0d0**0.5d0) do i=1,n so(i,1)=0.0d0 so(i,2)=0.0d0 so(i,3)=1.0d0 enddo call normal(vnorm,n,so) call vector(n,so,vnorm) call output(so,n) c orthog c60 to so here! call orthog(vnorm,n,c60,so) call vector(n,c60,vnorm) call checkn(so,n) call dotprod(k0,n,so,c60) write(9,*)'dotp so c60',k0 call cv(so,n) write(12,*)'c60 dipole field after so' call output(c60,n) call cv(c60,n) call sumit(n,c60,sumo) c ********************************************************** c V10 (-xz, -yz, 1-zz) write(12,*)'V10' do i=1,n v10(i,1)=(-1.0d0*(x(i)*z(i))) v10(i,2)=(-1.0d0*(y(i)*z(i))) v10(i,3)=1.0d0-((z(i)*z(i))) enddo 30 format(f23.20) write(10,*)'V10 raw' call checko(n,v10,so) write(10,*)'v10 end raw' call output(v10,n) call normal(vnorm,n,v10) call vector(n,v10,vnorm) call orthog(vnorm,n,v10,so) call vector(n,v10,vnorm) call normal(vnorm,n,v10) call vector(n,v10,vnorm) c v10 orth to c60 here! call orthog(vnorm,n,c60,v10) call vector(n,c60,vnorm) call checkn(v10,n) call checko(n,v10,so) call dotprod(k1,n,v10,c60) write(9,*)'dotp v10 c60',k1 c call cv(v10,n) write(12,*)'c60 dipole field after v10' call output(c60,n) call cv(c60,n) call sumit(n,c60,sumo) c *********************************************************** c V20 3(-xz2, -yz2, z(1-z2) do i=1,n v20(i,1)=3.0d0*(-x(i)*(z(i)**2)) v20(i,2)=3.0d0*(-y(i)*(z(i)**2)) v20(i,3)=3.0d0*(z(i)*(1.0d0-(z(i)**2))) enddo write(12,*)'V20' c call output(v20,n) write(10,*)'V20 raw' call checko(n,v20,so) call checko(n,v20,v10) write(10,*)'end raw' call normal(vnorm,n,v20) call vector(n,v20,vnorm) call checkn(v20,n) call checko(n,v20,so) call checko(n,v20,v10) call dotprod(dotp,n,v20,c60) write(9,*)'dotp v20 c60',dotp c *********************************************************** c V30 do i=1,n v3k=(15.0d0*(z(i)**2.0d0))-3.0d0 v30(i,1)=-0.5d0*x(i)*z(i)*v3k v30(i,2)=-0.5d0*y(i)*z(i)*v3k v30(i,3)=0.5d0*(1.0d0-(z(i)**2.0d0))*v3k enddo write(12,*)'V30' call output(v30,n) write(10,*)'V30 raw' call dotprod(dotp,n,v30,c60) write(6,*)dotp call checko(n,v30,so) call checko(n,v30,v10) call checko(n,v30,v20) write(10,*)'v30 end raw' call orthog(vnorm,n,v30,so) call vector(n,v30,vnorm) call normal(vnorm,n,v30) call vector(n,v30,vnorm) call orthog(vnorm,n,v30,v10) call vector(n,v30,vnorm) call normal(vnorm,n,v30) call vector(n,v30,vnorm) c make orthog v30 to c60 here! call orthog(vorth,n,c60,v30) call vector(n,c60,vorth) call checko(n,v30,so) call checko(n,v30,v10) call checko(n,v30,v20) call dotprod(k3,n,v30,c60) write(9,*)'dotp v30 c60',k3 c call cv(v30,n) write(12,*)'c60 dipole field after v30' call output(c60,n) call cv(c60,n) call sumit(n,c60,sumo) c *********************************************************** c V40 v4k=0.0d0 do i=1,n v4k=((35.0d0*(z(i)**3.0d0))-(15.0d0*z(i))) v40(i,1)=-0.5d0*x(i)*z(i)*v4k v40(i,2)=-0.5d0*y(i)*z(i)*v4k v40(i,3)=0.5d0*(1.0d0-(z(i)*z(i)))*v4k enddo c write(12,*)'V40' c call output(v40,n) write(10,*)'v40 raw' call checko(n,v40,so) call checko(n,v40,v10) call checko(n,v40,v20) call checko(n,v40,v30) write(10,*)'v40 raw end' call normal(vnorm,n,v40) call vector(n,v40,vnorm) c call orthog(vnorm,n,c60,v40) c call vector(n,c60,vnorm) c call checkn(c60,n) c call checko(n,c60,v40) c call orthog(vnorm,n,v40,v20) c call vector(n,v40,vnorm) c call checkn(v40,n) call checko(n,v40,so) call checko(n,v40,v10) call checko(n,v40,v20) call checko(n,v40,v30) call dotprod(dotp,n,v40,c60) write(9,*)'dotp v40 c60',dotp c *********************************************************** c V50 do i=1,n v5k=(315.0d0*(z(i)**4.0d0))-(210.0d0*(z(i)**2.0d0))+15.0d0 v50(i,1)=-0.125d0*x(i)*z(i)*v5k v50(i,2)=-0.125d0*y(i)*z(i)*v5k v50(i,3)=0.125d0*(1.0d0-(z(i)*z(i)))*v5k enddo write(10,*)'V50raw' write(12,*)'V50' call output(v50,n) call checko(n,v50,so) call checko(n,v50,v10) call checko(n,v50,v20) call checko(n,v50,v30) call checko(n,v50,v40) write(10,*)'v50 end raw' call orthog(vnorm,n,v50,so) call vector(n,v50,vnorm) call normal(vnorm,n,v50) call vector(n,v50,vnorm) call orthog(vnorm,n,v50,v10) call vector(n,v50,vnorm) call normal(vnorm,n,v50) call vector(n,v50,vnorm) call orthog(vnorm,n,v50,v30) call vector(n,v50,vnorm) call normal(vnorm,n,v50) call vector(n,v50,vnorm) c make orthog v50 to c60 here! call orthog(vorth,n,c60,v50) call vector(n,c60,vorth) call checko(n,v50,so) call checko(n,v50,v10) call checko(n,v50,v20) call checko(n,v50,v30) call checko(n,v50,v40) call dotprod(k5,n,v50,c60) write(9,*)'dotp v50 c60',k5 c call cv(v50,n) write(12,*)'c60 dipole field after v50' call output(c60,n) c c60(24,3)=0.0d0 c c60(25,3)=0.0d0 c c60(29,3)=0.0d0 c c60(32,3)=0.0d0 c c60(35,3)=0.0d0 c c60(36,3)=0.0d0 c c60(37,3)=0.0d0 c c60(56,3)=0.0d0 c write(12,*)'c60 dipole field after v50' c call output(c60,n) call cv(c60,n) call sumit(n,c60,sumo) c ******************************************************* c v60 write(12,*)'V60' do i=1,n v6k=((693.0d0*(z(i)**5.0d0))-(630.0d0*(z(i)**3.0d0))) v6k=v6k+(105.0d0*z(i)) v60(i,1)=(-0.125d0*x(i)*z(i))*v6k v60(i,2)=(-0.125d0*y(i)*z(i))*v6k v60(i,3)=(0.125d0*(1.0d0-(z(i)**2)))*v6k enddo c call output(v60,n) write(10,*)'V60 raw' call checko(n,v60,so) call checko(n,v60,v10) call checko(n,v60,v20) call checko(n,v60,v30) call checko(n,v60,v40) call checko(n,v60,v50) write(10,*)'v60 end raw' call orthog(vnorm,n,v60,v20) call vector(n,v60,vnorm) call checkn(v60,n) call orthog(vnorm,n,v60,v40) call vector(n,v60,vnorm) call checkn(v60,n) call checko(n,v60,so) call checko(n,v60,v10) call checko(n,v60,v20) call checko(n,v60,v30) call checko(n,v60,v40) call checko(n,v60,v50) call dotprod(dotp,n,v60,c60) write(9,*)'dotp v60 c60',dotp c *********************************************************** c v70 write(12,*)'V70' do i=1,n v7k=(3003.0d0*(z(i)**6.0d0))-(3465.0d0*(z(i)**4.0d0)) v7k=v7k+(945.0d0*(z(i)**2.0d0))-35.0d0 v70(i,1)=-0.0625d0*x(i)*z(i)*v7k v70(i,2)=-0.0625d0*y(i)*z(i)*v7k v70(i,3)=(0.0625d0*(1.0d0-z(i)**2))*v7k enddo call output(v70,n) write(10,*)'V70 raw' call checko(n,v70,so) call checko(n,v70,v10) call checko(n,v70,v20) call checko(n,v70,v30) call checko(n,v70,v40) call checko(n,v70,v50) call checko(n,v70,v60) write(10,*)'v70 end raw' call orthog(vnorm,n,v70,so) call vector(n,v70,vnorm) call normal(vnorm,n,v70) call vector(n,v70,vnorm) call orthog(vnorm,n,v70,v10) call vector(n,v70,vnorm) call normal(vnorm,n,v70) call vector(n,v70,vnorm) call orthog(vnorm,n,v70,v30) call vector(n,v70,vnorm) call normal(vnorm,n,v70) call vector(n,v70,vnorm) call orthog(vnorm,n,v70,v50) call vector(n,v70,vnorm) call normal(vnorm,n,v70) call vector(n,v70,vnorm) c set to orthog c60 to v70 call orthog(vnorm,n,c60,v70) call vector(n,c60,vnorm) call checko(n,v70,so) call checko(n,v70,v10) call checko(n,v70,v20) call checko(n,v70,v30) call checko(n,v70,v40) call checko(n,v70,v50) call checko(n,v70,v60) call dotprod(k7,n,v70,c60) write(9,*)'dotp v70 c60',k7 c call cv(v70,n) write(12,*)'c60 dipole field after v70' call output(c60,n) call cv(c60,n) call sumit(n,c60,sumo) c *********************************************************** c v80 write(12,*)'V80' do i=1,n v8k=(6435*(z(i)**7))-(9009*(z(i)**5))+(3465*(z(i)**3)) v8k=v8k-(315*z(i)) v80(i,1)=-0.0625d0*x(i)*z(i)*v8k v80(i,2)=-0.0625d0*y(i)*z(i)*v8k v80(i,3)=(0.0625d0*(1.0d0-z(i)**2))*v8k enddo call output(v80,n) write(10,*)'V80 raw' call checko(n,v80,so) call checko(n,v80,v10) call checko(n,v80,v20) call checko(n,v80,v30) call checko(n,v80,v40) call checko(n,v80,v50) call checko(n,v80,v60) call checko(n,v80,v70) write(10,*)'v80 end raw' call orthog(vnorm,n,v80,v20) call vector(n,v80,vnorm) call checkn(v80,n) call orthog(vnorm,n,v80,v40) call vector(n,v80,vnorm) call checkn(v80,n) call orthog(vnorm,n,v80,v60) call vector(n,v80,vnorm) call checkn(v80,n) call checko(n,v80,so) call checko(n,v80,v10) call checko(n,v80,v20) call checko(n,v80,v30) call checko(n,v80,v40) call checko(n,v80,v50) call checko(n,v80,v60) call checko(n,v80,v70) call dotprod(dotp,n,v80,c60) write(9,*)'dotp v80 c60',dotp c ************************************************************ c v90 write(12,*)'v90' c x(1)=4.0d0 c y(1)=3.0d0 c z(1)=2.0d0 do i=1,n v9k=(109395.0d0*(z(i)**8.0d0))-(180180.0d0*(z(i)**6.0d0)) v9k=v9k+(90090.0d0*(z(i)**4.0d0)) v9k=v9k-(13860.0d0*(z(i)**2.0d0))+315.0d0 v90(i,1)=-(1.0d0/128.0d0)*x(i)*z(i)*v9k v90(i,2)=-(1.0d0/128.0d0)*y(i)*z(i)*v9k v90(i,3)=(1.0d0/128.0d0)*(1.0d0-z(i)**2)*v9k enddo c write(6,*)v90(1,1),v90(1,2),v90(1,3) call output(v90,n) write(10,*)'v90 raw' call checko(n,v90,so) call checko(n,v90,v10) call checko(n,v90,v20) call checko(n,v90,v30) call checko(n,v90,v40) call checko(n,v90,v50) call checko(n,v90,v60) call checko(n,v90,v70) call checko(n,v90,v80) write(10,*)'v90 end raw' call orthog(vnorm,n,v90,so) call vector(n,v90,vnorm) call normal(vnorm,n,v90) call vector(n,v90,vnorm) call orthog(vnorm,n,v90,v10) call vector(n,v90,vnorm) call normal(vnorm,n,v90) call vector(n,v90,vnorm) call orthog(vnorm,n,v90,v30) call vector(n,v90,vnorm) call normal(vnorm,n,v90) call vector(n,v90,vnorm) call orthog(vnorm,n,v90,v50) call vector(n,v90,vnorm) call normal(vnorm,n,v90) call vector(n,v90,vnorm) call orthog(vnorm,n,v90,v70) call vector(n,v90,vnorm) call normal(vnorm,n,v90) call vector(n,v90,vnorm) c orthog v90 to c60 call orthog(vnorm,n,c60,v90) call vector(n,c60,vnorm) call checko(n,v90,so) call checko(n,v90,v10) call checko(n,v90,v20) call checko(n,v90,v30) call checko(n,v90,v40) call checko(n,v90,v50) call checko(n,v90,v60) call checko(n,v90,v70) call checko(n,v90,v80) call dotprod(k9,n,v90,c60) write(9,*)'dotp v90 c60',k9 c call cv(v90,n) write(12,*)'c60 dipole field after v90' call output(c60,n) call cv(c60,n) call sumit(n,c60,sumo) c ************************************************************ c v100 write(12,*)'v100' do i=1,n v100k=(230945*z(i)**9)-(437580*z(i)**7)+(270270*z(i)**5) v100k=v100k-(60060*(z(i)**3))+(3465*(z(i))) v100(i,1)=-(1.0d0/128.0d0)*x(i)*z(i)*v100k v100(i,2)=-(1.0d0/128.0d0)*y(i)*z(i)*v100k v100(i,3)=(1.0d0/128.0d0)*(1.0d0-z(i)**2)*v100k enddo call output(v100,n) write(10,*)'v100 raw' call checko(n,v100,so) call checko(n,v100,v10) call checko(n,v100,v20) call checko(n,v100,v30) call checko(n,v100,v40) call checko(n,v100,v50) call checko(n,v100,v60) call checko(n,v100,v70) call checko(n,v100,v80) call checko(n,v100,v90) write(10,*)'end raw' call dotprod(dotp,n,v100,c60) write(9,*)'dotp v100 c60',dotp c ************************************************************ c v110 write(12,*)'v110' do i=1,n v110k=(969969.0d0*(z(i)**10.0d0))-(2078505.0d0*(z(i)**8)) v110k=v110k+(1531530.0d0*(z(i)**6))-(450450.0d0*(z(i)**4)) v110k=v110k+(45045.0d0*(z(i)**2.0d0))-693.0d0 v110(i,1)=(-1.0d0/256.0d0)*x(i)*z(i)*v110k v110(i,2)=(-1.0d0/256.0d0)*y(i)*z(i)*v110k v110(i,3)=(1.0d0/256.0d0)*(1.0d0-z(i)**2)*v110k enddo call output(v110,n) write(12,*)'c60 upto 110' call output(c60,n) write(10,*)'v110 raw' call checko(n,v110,so) call checko(n,v110,v10) call checko(n,v110,v20) call checko(n,v110,v30) call checko(n,v110,v40) call checko(n,v110,v50) call checko(n,v110,v60) call checko(n,v110,v70) call checko(n,v110,v80) call checko(n,v110,v90) call checko(n,v110,v100) write(10,*)'end raw' call orthog(vnorm,n,v110,so) call vector(n,v110,vnorm) call normal(vnorm,n,v110) call vector(n,v110,vnorm) call orthog(vnorm,n,v110,v10) call vector(n,v110,vnorm) call normal(vnorm,n,v110) call vector(n,v110,vnorm) call orthog(vnorm,n,v110,v30) call vector(n,v110,vnorm) call normal(vnorm,n,v110) call vector(n,v110,vnorm) call orthog(vnorm,n,v110,v50) call vector(n,v110,vnorm) call normal(vnorm,n,v110) call vector(n,v110,vnorm) call orthog(vnorm,n,v110,v70) call vector(n,v110,vnorm) call normal(vnorm,n,v110) call vector(n,v110,vnorm) call orthog(vnorm,n,v110,v90) call vector(n,v110,vnorm) call normal(vnorm,n,v110) call vector(n,v110,vnorm) c make orthog here call orthog(vnorm,n,c60,v110) call vector(n,c60,vnorm) call checko(n,v110,so) call checko(n,v110,v10) call checko(n,v110,v20) call checko(n,v110,v30) call checko(n,v110,v40) call checko(n,v110,v50) call checko(n,v110,v60) call checko(n,v110,v70) call checko(n,v110,v80) call checko(n,v110,v90) call checko(n,v110,v100) call dotprod(k11,n,v110,c60) write(9,*)'dotp v110 c60',k11 c call cv(c60,n) write(12,*)'c60 dipole field after v90' call output(c60,n) call cv(c60,n) call sumit(n,c60,sumo) c *********************************************************** c v130 write(12,*)'v130' do i=1,n v130k=(16900975*(z(i)**12))-(44618574*(z(i)**10)) v130k=v130k+(43648605*(z(i)**8))-(19399380*(z(i)**6)) v130k=v130k+(3828825*(z(i)**4))-(270270*(z(i)**2)) v130k=v130k+3003 v130(i,1)=-(1.0d0/1024.0d0)*x(i)*z(i)*v130k v130(i,2)=-(1.0d0/1024.0d0)*y(i)*z(i)*v130k v130(i,3)=(1.0d0/1024.0d0)*(1.0d0-z(i)**2)*v130k enddo call output(v130,n) write(10,*)'v130 raw' call checko(n,v130,so) call checko(n,v130,v10) call checko(n,v130,v20) call checko(n,v130,v30) call checko(n,v130,v40) call checko(n,v130,v50) call checko(n,v130,v60) call checko(n,v130,v70) call checko(n,v130,v80) call checko(n,v130,v90) call checko(n,v130,v100) call checko(n,v130,v110) write(10,*)'end raw' c make orthog here c call orthog(vnorm,n,c60,v130) c call vector(n,c60,vnorm) call orthog(vnorm,n,v130,so) call vector(n,v130,vnorm) call normal(vnorm,n,v130) call vector(n,v130,vnorm) call orthog(vnorm,n,v130,v10) call vector(n,v130,vnorm) call normal(vnorm,n,v130) call vector(n,v130,vnorm) call orthog(vnorm,n,v130,v30) call vector(n,v130,vnorm) call normal(vnorm,n,v130) call vector(n,v130,vnorm) call orthog(vnorm,n,v130,v50) call vector(n,v130,vnorm) call normal(vnorm,n,v130) call vector(n,v130,vnorm) call orthog(vnorm,n,v130,v70) call vector(n,v130,vnorm) call normal(vnorm,n,v130) call vector(n,v130,vnorm) call orthog(vnorm,n,v130,v90) call vector(n,v130,vnorm) call normal(vnorm,n,v130) call vector(n,v130,vnorm) call orthog(vnorm,n,v130,v110) call vector(n,v130,vnorm) call normal(vnorm,n,v130) call vector(n,v130,vnorm) call checko(n,v130,so) call checko(n,v130,v10) call checko(n,v130,v20) call checko(n,v130,v30) call checko(n,v130,v40) call checko(n,v130,v50) call checko(n,v130,v60) call checko(n,v130,v70) call checko(n,v130,v80) call checko(n,v130,v90) call checko(n,v130,v100) call checko(n,v130,v110) call dotprod(k13,n,v130,c60) write(9,*)'dotp v130 c60',k13 call cv(v130,n) c ************************************************************ c v150 here write(12,*)'v150' do i=1,n v150k=(145422675*(z(i)**14))-(456326325*(z(i)**12)) v150k=v150k+(557732175*(z(i)**10))-(334639305*(z(i)**8)) v150k=v150k+(101846745*(z(i)**6))-(14549535*(z(i)**4)) v150k=v150k+(765765*(z(i)**2))-(6435) v150(i,1)=-(1.0d0/2048.0d0)*x(i)*z(i)*v150k v150(i,2)=-(1.0d0/2048.0d0)*y(i)*z(i)*v150k v150(i,3)=(1.0d0/2048.0d0)*(1.0d0-z(i)**2)*v150k enddo call output(v150,n) call normal(vnorm,n,v150) call vector(n,v150,vnorm) write(10,*)'v150 raw' call checko(n,v150,so) call checko(n,v150,v10) c call checko(n,v150,v20) call checko(n,v150,v30) c call checko(n,v150,v40) call checko(n,v150,v50) c call checko(n,v150,v60) call checko(n,v150,v70) c call checko(n,v150,v80) call checko(n,v150,v90) c call checko(n,v150,v100) call checko(n,v150,v110) call checko(n,v150,v130) write(10,*)'end raw' call normal(vnorm,n,v150) call vector(n,v150,vnorm) call orthog(vnorm,n,v150,so) call vector(n,v150,vnorm) call normal(vnorm,n,v150) call vector(n,v150,vnorm) call orthog(vnorm,n,v150,v10) call vector(n,v150,vnorm) call normal(vnorm,n,v150) call vector(n,v150,vnorm) call orthog(vnorm,n,v150,v20) call vector(n,v150,vnorm) call normal(vnorm,n,v150) call vector(n,v150,vnorm) call orthog(vnorm,n,v150,v30) call vector(n,v150,vnorm) call normal(vnorm,n,v150) call vector(n,v150,vnorm) call orthog(vnorm,n,v150,v50) call vector(n,v150,vnorm) call normal(vnorm,n,v150) call vector(n,v150,vnorm) call orthog(vnorm,n,v150,v70) call vector(n,v150,vnorm) call normal(vnorm,n,v150) call vector(n,v150,vnorm) call orthog(vnorm,n,v150,v90) call vector(n,v150,vnorm) call normal(vnorm,n,v150) call vector(n,v150,vnorm) call orthog(vnorm,n,v150,v110) call vector(n,v150,vnorm) call normal(vnorm,n,v150) call vector(n,v150,vnorm) call orthog(vnorm,n,v150,v130) call vector(n,v150,vnorm) call normal(vnorm,n,v150) call vector(n,v150,vnorm) call checko(n,v150,so) call checko(n,v150,v10) c call checko(n,v150,v20) call checko(n,v150,v30) c call checko(n,v150,v40) call checko(n,v150,v50) c call checko(n,v150,v60) call checko(n,v150,v70) c call checko(n,v150,v80) call checko(n,v150,v90) c call checko(n,v150,v100) call checko(n,v150,v110) call checko(n,v150,v130) call dotprod(k15,n,v150,c60) write(9,*)'dotp v150 c60',k15 call cv(v150,n) c ************************************************************ c Calculations c a(cage)=coeff(scalar).vsh sum off etc do i=1,n do j=1,3 cage(i,j)=(k0*so(i,j))+(k1*v10(i,j)) cage(i,j)=cage(i,j)+(k5*v50(i,j))+(k7*v70(i,j)) cage(i,j)=cage(i,j)+(k9*v90(i,j)) c cage(i,j)=cage(i,j)+(k11*v110(i,j)) c cage(i,j)=cage(i,j)+(k13*v130(i,j)) cage(i,j)=cage(i,j)+(k15*v150(i,j)) c cage(i,j)=cage(i,j)+(k17*v170(i,j)) enddo enddo c 25 format(4x,d20.12,1x,d20.12,1x,d20.12) do i=1,n x(i)=cage(i,1) y(i)=cage(i,2) z(i)=cage(i,3) write(14,25)x(i),y(i),z(i) enddo call orthog(vnorm,n,c60,cage) call vector(n,c60,vnorm) c call normal(vnorm,n,c60) c call vector(n,c60,vnorm) write(12,*)'c60 orthog cage' call output(c60,n) c testing how far away as same vector call dotprod(dotp,n,c60,c60) write(14,*)'dot c60 c60',dotp call dotprod(dotp,n,c60,cage) write(14,*)'dot c60 cage',dotp c ************************************************************ stop end subroutine dotprod(dotp, n, a, b) implicit real*8(a-h,o-z) integer n double precision dotp dimension a(n,3), b(n,3) dotp=0.0d0 do i=1,n do j=1,3 dotp=dotp+(a(i,j)*b(i,j)) enddo enddo return end c __________________________________________________________ subroutine normal(vnorm,n,a) implicit real*8(a-h,o-z) integer n double precision dotp,sqrk dimension vnorm(n,3),a(n,3) call dotprod(dotp,n,a,a) do i=1,n do j=1,3 sqrk=SQRT(abs(dotp)) vnorm(i,j)=a(i,j)/sqrk enddo enddo return end c __________________________________________________________ subroutine orthog(c,n,a,b) implicit real*8(a-h,o-z) integer n double precision dotp dimension vorth(n,3),a(n,3),b(n,3),c(n,3) call dotprod(dotp,n,a,b) do i=1,n do j=1,3 vorth(i,j)=a(i,j)-(dotp*b(i,j)) enddo enddo c call normal(c,n,vorth) call vector(n,c,vorth) return end c __________________________________________________________ subroutine output(vout, n) implicit real*8(a-h,o-z) integer n dimension vout(n,3) 99 format(I2,1x,f24.19,1x,f24.19,1x,f24.19) do i=1,n write(12,99)i,vout(i,1), vout(i,2), vout(i,3) enddo return end c __________________________________________________________ subroutine checkn(vnorm, n) implicit real*8(a-h,o-z) integer n double precision dotp dimension vnorm(n,3) call dotprod(dotp, n, vnorm, vnorm) write(10,*)'Normalisation this should be = 1',dotp return end c _________________________________________________________ subroutine checko(n,vx,vy) implicit real*8(a-h,o-z) integer n double precision dotp dimension vx(n,3),vy(n,3) call dotprod(dotp,n,vx,vy) write(10,*)'Orthogalisation this should be = 0',dotp return end c _________________________________________________________ subroutine vector(n,a,b) implicit real*8(a-h,o-z) integer n dimension a(n,3),b(n,3) do i=1,n do j=1,3 a(i,j)=b(i,j) enddo enddo return end c _________________________________________________________ subroutine cv(x,n) implicit real*8(a-h,o-z) integer n dimension x(n,3),sum(3) do i=1,n sum(1)=sum(1)+x(i,1) sum(2)=sum(2)+x(i,2) sum(3)=sum(3)+x(i,3) enddo write(12,*)'x',sum(1),' y',sum(2),' z',sum(3) return end c ___________________________________________________ subroutine sumit(n,a,sumo) implicit real*8(a-h,o-z) integer n double precision sumo dimension a(n,3) sumo=0.0d0 do i=1,n do j=1,3 sumo=sumo+abs(a(i,j)) enddo enddo write(12,*)'sumo',sumo return end