PROGRAM PESWF
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
parameter (nvar0=4,nh0=4100,nbas0=100,ngr0=400,nhdim0=4100,
1ngp0=140000)
logical lper(nvar0),lprop,lh,ltime,lwr
character*5 s5
dimension ngr(nvar0),ngrh(nvar0),nbase(nvar0),nh(nvar0),
1nbas(nvar0),ic(nvar0),ice(nvar0),ich(nvar0),l1(ngp0),l2(ngp0),
3l3(ngp0),ip(nh0),peri(nvar0)
real*8 min(nvar0),max(nvar0),G11k(nbas0),enk(nbas0),
1ens(ngp0),G11s(ngr0),xs(ngr0),enkl(nbas0,nbas0),pns(ngp0),
2G11kl(nbas0,nbas0),G12kl(nbas0,nbas0),G22kl(nbas0,nbas0),
3Gkl(nbas0,nbas0),var(ngp0,nvar0),kelvin,dre,anorm
common//h(nhdim0,nhdim0),c(nhdim0,nhdim0),q(nhdim0),
1s(nhdim0,nhdim0),vu(nhdim0,nhdim0),d(nhdim0),
1dn(nhdim0,nhdim0),popul(nhdim0),pa(nhdim0),pa2(nhdim0),
1dre(nhdim0)
common/const/pi,bohr
common/ho/p(nh0,nh0),ank(nvar0,nh0),w(nvar0),x0(nvar0),am0(nvar0)
common/pw/kindpw
real*8 dt,xt1,xt2,htocm
integer*4 ntstep
common/timec/dt,xt1,xt2,ntstep
c
pi=4.0d0*atan(1.0d0)
bohr=0.529177d0
c hartree to cm-1:
htocm=219474.63d0
c
c angles supposed in deg (o)
c other coords in A
c
write(6,6000)
6000 format(/,' Molecular PESs',/,/,
1' This program reads PES and metric tensor in non-orthogonal',/,
1' coordinates and calculates energies',/,/,
1 ' Input: SURF.OPT ... options',/,
1 ' ENN.LST ... list of variables and energies',/,
1 ' PROPERTYN.LST ... list of properties (Opt)',/)
C
call getpar(nvar,nvar0,lper,ngr,min,max,nbase,nbas,nh,ngrh,
1ice,ic,ich,wcm,x0,am0,ipl,ip,ikin,kindpw,kinddia,kelvin,
2w,nbas0,nh0,lh,emax1,emax2,peri,ltime,lwr)
call check(ipl,ikin,nvar,ic,kindpw)
c
io=0
do 24 ii=1,nvar
24 if(ic(ii).eq.1.or.ich(ii).eq.1.or.ice(ii).eq.1)io=1
if(io.eq.1)then
maxb=1
do 25 ii=1,nvar
if(nbas(ii).gt.maxb.and.ic(ii).eq.1)maxb=nbas(ii)
if(nh(ii).gt.maxb.and.ich(ii).eq.1)maxb=nh(ii)
25 if(nbase(ii).gt.maxb.and.ice(ii).eq.1)maxb=nbase(ii)
write(6,*)'Initializing HO basis, kmax = ',maxb
do 103 ii=1,maxb
do 103 jj=1,maxb
103 p(jj,ii)=0.0d0
p(1,1)=1.0d0
do 105 ii=1,nvar
105 ank(ii,1)=dsqrt(dsqrt(am0(ii)*w(ii)/pi)/bohr)
k=1
102 k=k+1
akfac=1.0d0
do 101 kk=1,k-1
101 akfac=akfac*dble(kk)
do ii=1,nvar
ank(ii,k)=dsqrt(dsqrt(am0(ii)*w(ii)/pi)/bohr/2.0d0**(k-1)/akfac)
enddo
p(k,1 )=-p(k-1,2)
do 104 j=2,k-2
104 p(k,j )=(2.0d0*p(k-1,j-1)-dble(j)*p(k-1,j+1))
if(k-2.gt.0)p(k,k-1)=2.0d0*p(k-1,k-2)
p(k,k )=2.0d0*p(k-1,k-1)
if(k.lt.maxb)goto 102
endif
c
do 11 k=1,nbase(1)
enk(k)=0.0d0
do 11 l=1,nbase(2)
11 enkl(k,l)=0.0d0
do 1 k=1,nbas(1)
G11k(k)=0.0d0
do 1 l=1,nbas(2)
Gkl(k,l)=0.0d0
G22kl(k,l)=0.0d0
G12kl(k,l)=0.0d0
1 G11kl(k,l)=0.0d0
c read grid and calculate expansion in the basis
open(8,file='ENN.LST')
inquire(file='PROPERTYN.LST',exist=lprop)
if(lprop)open(86,file='PROPERTYN.LST')
write(6,*)'Reading E,G'
if(ngr(1)+1.gt.ngr0)call report('Too many grid points')
al=max(1)-min(1)
dx=al/dble(ngr(1))
x=min(1)-dx
if(nvar.eq.1)then
ie=0
do 2 ix=0,ngr(1)
x=x+dx
ie=ie+1
read(8,*)idumm,xdumm,en,G11
if(lprop)read(86,*)xdumm,propn
if(lprop)pns(ix+1)=propn
xs(ix+1)=x
ens(ix+1)=en
G11s(ix+1)=G11
if(ix.eq.0.or.ix.eq.ngr(1))then
dd=dx/2.0d0
else
dd=dx
endif
do 22 k=1,nbase(1)
22 enk(k) =enk(k) +en *fik(1,ic(1),0,k,x,min(1),al)*dd
do 2 k=1,nbas(1)
2 G11k(k)=G11k(k)+G11*fik(1,ic(1),0,k,x,min(1),al)*dd
else
if(nvar.eq.2)then
bl=max(2)-min(2)
dy=bl/dble(ngr(2))
ie=0
do 3 ix=0,ngr(1)
x=x+dx
y=min(2)-dy
if(ix.eq.0.or.ix.eq.ngr(1))then
ddx=dx/2.0d0
else
ddx=dx
endif
do 3 iy=0,ngr(2)
ie=ie+1
y=y+dy
if(iy.eq.0.or.iy.eq.ngr(2))then
ddy=dy/2.0d0
else
ddy=dy
endif
ds=ddy*ddx
read(8,*)idum,idum,xdumm,ydumm,en,G,G11,G12,G22
if(ie.gt.ngp0)call report('Too many grid points')
ens(ie)=en
if(lprop)read(86,*)xdumm,ydumm,propn
if(lprop)pns(ie)=propn
var(ie,1)=x
var(ie,2)=y
do 422 k=1,nbase(1)
fk=fik(1,ic(1),0,k,x,min(1),al)*en*ds
do 422 l=1,nbase(2)
422 enkl(k,l) =enkl(k,l) +fk*fik(2,ic(2),0,l,y,min(2),bl)
do 3 k=1,nbas(1)
fk=fik(1,ic(1),0,k,x,min(1),al)*ds
do 3 l=1,nbas(2)
fkl=fk*fik(2,ic(2),0,l,y,min(2),bl)
Gkl(k,l)=Gkl(k,l)+G*fkl
G22kl(k,l)=G22kl(k,l)+G22*fkl
G12kl(k,l)=G12kl(k,l)+G12*fkl
3 G11kl(k,l)=G11kl(k,l)+G11*fkl
else
call report('higher dimensions not implemented')
endif
endif
close(8)
if(lprop)close(86)
c
c Reproduce ENN:
open(8,file='ENNT.LST')
x=min(1)-dx
if(nvar.eq.1)then
do 21 ix=0,ngr(1)
x=x+dx
if(ice(1).ne.1.and.ice(1).ne.0)then
en=ens(ix+1)
else
en=0.0d0
do 222 k=1,nbase(1)
222 en =en +enk(k)*fik(1,ic(1),0,k,x,min(1),al)
endif
if(ic(1).ne.1.and.ic(1).ne.0)then
G11=G11s(ix+1)
else
G11=0.0d0
do 221 k=1,nbas(1)
221 G11=G11+G11k(k)*fik(1,ic(1),0,k,x,min(1),al)
endif
21 write(8,808)x,en,G11
endif
if(nvar.eq.2)then
ie=0
do 31 ix=0,ngr(1)
x=x+dx
y=min(2)-dy
do 31 iy=0,ngr(2)
ie=ie+1
y=y+dy
if(ice(1).ne.1.and.ice(1).ne.0)then
en=ens(ie)
else
en=0.0d0
do 322 k=1,nbase(1)
fk=fik(1,ic(1),0,k,x,min(1),al)
do 322 l=1,nbase(2)
322 en=en+enkl(k,l)*fk*fik(2,ic(2),0,l,y,min(2),bl)
endif
G11=0.0d0
G12=0.0d0
G22=0.0d0
G=0.0d0
do 321 k=1,nbas(1)
fk=fik(1,ic(1),0,k,x,min(1),al)
do 321 l=1,nbas(2)
fkl=fk*fik(2,ic(2),0,l,y,min(2),bl)
G=G+Gkl(k,l)*fkl
G11=G11+G11kl(k,l)*fkl
G12=G12+G12kl(k,l)*fkl
321 G22=G22+G22kl(k,l)*fkl
31 write(8,809)x,y,en,G,G11,G12,G22
809 format(3f12.6,4g13.5)
endif
close(8)
nhdim=1
do 27 ii=1,nvar
27 nhdim=nhdim*nh(ii)
write(6,6004)nhdim
6004 format(' Dimension of the Hamiltonian = ',i6)
if(nhdim.gt.nhdim0)call report(' ... too high.')
do 7 i=1,nhdim
do 7 j=1,nhdim
s(i,j)=0.0d0
7 h(i,j)=0.0d0
t22=0.0d0
v22=0.0d0
t11=0.0d0
v11=0.0d0
c
write(6,*)'Making H'
open(81,file='ENNTH.LST')
if(lprop)open(816,file='PROPERTYNTH.LST')
if(nvar.eq.1)then
dxh=al/dble(ngrh(1))
x=min(1)-dxh
do 217 ix=0,ngrh(1)
if(ix.eq.0.or.ix.eq.ngrh(1))then
dd=dxh/2.0d0
else
dd=dxh
endif
x=x+dxh
if(ic(1).eq.0.or.ic(1).eq.1)then
G=0.0d0
Gq=0.0d0
Gqq=0.0d0
do 2217 k=1,nbas(1)
G =G +G11k(k)*fik(1,ic(1),0,k,x,min(1),al)
Gq =Gq +G11k(k)*fik(1,ic(1),1,k,x,min(1),al)
2217 Gqq=Gqq+G11k(k)*fik(1,ic(1),2,k,x,min(1),al)
else
call find(i1,i2,i3,xs,x,lper(1),ngr(1)+1,peri(1))
call abc(a,b,cc,xs(i1),xs(i2),xs(i3),G11s(i1),G11s(i2),G11s(i3))
G=x*(a*x+b)+cc
Gq=2.0d0*a*x+b
Gqq=2.0d0*a
endif
if(ice(1).eq.0.or.ice(1).eq.1)then
en=0.0d0
do 22171 k=1,nbase(1)
22171 en =en +enk(k) *fik(1,ic(1),0,k,x,min(1),al)
else
call find(i1,i2,i3,xs,x,lper(1),ngr(1)+1,peri(1))
call abc(a,b,cc,xs(i1),xs(i2),xs(i3),ens(i1),ens(i2),ens(i3))
en=x*(a*x+b)+cc
endif
write(81,808)x,en,G,Gq,Gqq
808 format(2f12.6,5g13.5)
if(lprop)then
call find(i1,i2,i3,xs,x,lper(1),ngr(1)+1,peri(1))
call abc(a,b,cc,xs(i1),xs(i2),xs(i3),pns(i1),pns(i2),pns(i3))
pn=x*(a*x+b)+cc
write(816,808)x,pn
endif
c
do 217 i=1,nh(1)
fi =fik(1,ich(1),0,i,x,min(1),al)
fiq =fik(1,ich(1),1,i,x,min(1),al)
do 217 j=1,nh(1)
fj =fik(1,ich(1),0,j,x,min(1),al)
fjq =fik(1,ich(1),1,j,x,min(1),al)
fjqq=fik(1,ich(1),2,j,x,min(1),al)
if(dabs(G).gt.1.0d-30)then
if(ikin.eq.1)then
t=0.125d0*G*(4.0d0*fiq*fjq+0.25d0/G**2*Gq**2*fi*fj
1 +(fiq*fj+fi*fjq)*Gq/G)*dd
else
t=-fi*0.125d0*((-0.25d0*Gq**2/G+Gqq)*fj
1 +4.0d0*Gq*fjq+4.0d0*G*fjqq)*dd
endif
endif
v=fi*en*fj*dd
if(i.eq.2.and.j.eq.2)then
t22=t22+t
v22=v22+v
endif
if(i.eq.1.and.j.eq.1)then
t11=t11+t
v11=v11+v
endif
s(i,j)=s(i,j)+fi*fj*dd
217 h(i,j)=h(i,j)+t+v
endif
if(nvar.eq.2)then
dxh=al/dble(ngrh(1))
dyh=bl/dble(ngrh(2))
if(lh)then
write(6,5600)emax1,emax2
5600 format(' Energy limit imposed: ',f12.6,'->',f12.6,' hartree')
do 756 i=1,nhdim
s(i,i)=1.0d0
756 h(i,i)=emax2
if(kinddia.eq.1)call report('Orthogonalization not allowed'//
1 ' with emax')
endif
x=min(1)-dxh
ipr=0
do 218 ix=0,ngrh(1)
if(ix*100/ngrh(1).gt.ipr)then
ipr=ix*100/ngrh(1)
write(6,6006)ipr
6006 format(i3,'%',$)
endif
x=x+dxh
if(ix.eq.0.or.ix.eq.ngrh(1))then
ddx=dxh/2.0d0
else
ddx=dxh
endif
y=min(2)-dyh
do 218 iy=0,ngrh(2)
y=y+dyh
if(iy.eq.0.or.iy.eq.ngrh(2))then
ddy=dyh/2.0d0
else
ddy=dyh
endif
ds=ddy*ddx
if(ic(1).eq.0.or.ic(1).eq.1)then
G=0.0d0
Gp=0.0d0
Gq=0.0d0
Gpq=0.0d0
G11=0.0d0
G12=0.0d0
G22=0.0d0
do 2218 k=1,nbas(1)
fk =fik(1,ic(1),0,k,x,min(1),al)
fkp =fik(1,ic(1),1,k,x,min(1),al)
do 2218 l=1,nbas(2)
fkl =fk *fik(2,ic(2),0,l,y,min(2),bl)
fklp =fkp*fik(2,ic(2),0,l,y,min(2),bl)
fklq =fk *fik(2,ic(2),1,l,y,min(2),bl)
fklpq=fkp*fik(2,ic(2),1,l,y,min(2),bl)
G =G +Gkl(k,l)*fkl
Gp =Gp +Gkl(k,l)*fklp
Gq =Gq +Gkl(k,l)*fklq
Gpq=Gpq+Gkl(k,l)*fklpq
G11=G11+G11kl(k,l)*fkl
G12=G12+G12kl(k,l)*fkl
2218 G22=G22+G22kl(k,l)*fkl
else
call report('Quadratic G-fit not implemented in 2D')
endif
if(ice(1).eq.0.or.ice(1).eq.1)then
if(lprop)call report('With property E must be with splines')
en=0.0d0
do 22181 k=1,nbase(1)
fk=fik(1,ic(1),0,k,x,min(1),al)
do 22181 l=1,nbase(2)
22181 en =en+enkl(k,l)*fk*fik(2,ic(2),0,l,y,min(2),bl)
else
call findlist(l1,l2,l3,var,x,lper(1),ie,ngp0,nvar0,peri(1))
q1=var(l1(2),1)
q2=var(l2(2),1)
q3=var(l3(2),1)
if(lper(1))then
tp=peri(1)
tph=tp/2.0d0
if(q1-x.gt.tph)q1=q1-tp
if(x-q1.gt.tph)q1=q1+tp
if(q2-x.gt.tph)q2=q2-tp
if(x-q2.gt.tph)q2=q2+tp
if(q3-x.gt.tph)q3=q3-tp
if(x-q3.gt.tph)q3=q3+tp
endif
c
c l1 l2 l3
c
c |
c |
c | p13
c p11 p12 |
c p22 |
c y-------------p21-----------+---------p23---- f1 f2 f3 values
c | from Ps for q1 q2 q3
c p32 |
c p31 | p33
c |
c |
c |
c |
c q1 q2 | q3 closest calculated points
c |
c |
c x x value
c
call findy(y,l1,i11,i21,i31,var,lper(2),ngp0,nvar0,peri(2))
p11=var(l1(i11),2)
p21=var(l1(i21),2)
p31=var(l1(i31),2)
call findy(y,l2,i12,i22,i32,var,lper(2),ngp0,nvar0,peri(2))
p12=var(l2(i12),2)
p22=var(l2(i22),2)
p32=var(l2(i32),2)
call findy(y,l3,i13,i23,i33,var,lper(2),ngp0,nvar0,peri(2))
p13=var(l3(i13),2)
p23=var(l3(i23),2)
p33=var(l3(i33),2)
if(lper(2))then
tp=peri(2)
tph=tp/2.0d0
if(p11-y.gt.tph)p11=p11-tp
if(y-p11.gt.tph)p11=p11+tp
if(p21-y.gt.tph)p21=p21-tp
if(y-p21.gt.tph)p21=p21+tp
if(p31-y.gt.tph)p31=p31-tp
if(y-p31.gt.tph)p31=p31+tp
if(p12-y.gt.tph)p12=p12-tp
if(y-p12.gt.tph)p12=p12+tp
if(p22-y.gt.tph)p22=p22-tp
if(y-p22.gt.tph)p22=p22+tp
if(p32-y.gt.tph)p32=p32-tp
if(y-p32.gt.tph)p32=p32+tp
if(p13-y.gt.tph)p13=p13-tp
if(y-p13.gt.tph)p13=p13+tp
if(p23-y.gt.tph)p23=p23-tp
if(y-p23.gt.tph)p23=p23+tp
if(p33-y.gt.tph)p33=p33-tp
if(y-p33.gt.tph)p33=p33+tp
endif
g1=ens(l1(i11))
g2=ens(l1(i21))
g3=ens(l1(i31))
call abc(a,b,cc,p11,p21,p31,g1,g2,g3)
f1=y*(a*y+b)+cc
g1=ens(l2(i12))
g2=ens(l2(i22))
g3=ens(l2(i32))
call abc(a,b,cc,p12,p22,p32,g1,g2,g3)
f2=y*(a*y+b)+cc
g1=ens(l3(i13))
g2=ens(l3(i23))
g3=ens(l3(i33))
call abc(a,b,cc,p13,p23,p33,g1,g2,g3)
f3=y*(a*y+b)+cc
call abc(a,b,cc,q1,q2,q3,f1,f2,f3)
en=x*(a*x+b)+cc
if(lprop)then
g1=pns(l1(i11))
g2=pns(l1(i21))
g3=pns(l1(i31))
call abc(a,b,cc,p11,p21,p31,g1,g2,g3)
f1=y*(a*y+b)+cc
g1=pns(l2(i12))
g2=pns(l2(i22))
g3=pns(l2(i32))
call abc(a,b,cc,p12,p22,p32,g1,g2,g3)
f2=y*(a*y+b)+cc
g1=pns(l3(i13))
g2=pns(l3(i23))
g3=pns(l3(i33))
call abc(a,b,cc,p13,p23,p33,g1,g2,g3)
f3=y*(a*y+b)+cc
call abc(a,b,cc,q1,q2,q3,f1,f2,f3)
pn=x*(a*x+b)+cc
write(816,809)x,y,pn
endif
endif
write(81,809)x,y,en,G,G11,G12,G22
qgpp=0.0d0
qgpq=0.0d0
qgqq=0.0d0
if(dabs(G).gt.0.0d0)then
qgpp=0.25d0/G**2*Gp*Gp
qgpq=0.25d0/G**2*Gp*Gq
qgqq=0.25d0/G**2*Gq*Gq
endif
c
if(.not.lh.or.en.lt.emax1)then
i=0
do 2181 ki=1,nh(1)
fk =fik(1,ich(1),0,ki,x,min(1),al)
fkp=fik(1,ich(1),1,ki,x,min(1),al)
do 2181 li=1,nh(2)
i=i+1
fl = fik(2,ich(2),0,li,y,min(2),bl)
fi =fk *fl
fip =fkp*fl
fiq =fk* fik(2,ich(2),1,li,y,min(2),bl)
j=0
do 2181 kj=1,nh(1)
fkj =fik(1,ich(1),0,kj,x,min(1),al)
fkjp=fik(1,ich(1),1,kj,x,min(1),al)
do 2181 lj=1,nh(2)
j=j+1
if(j.le.i)then
flj = fik(2,ich(2),0,lj,y,min(2),bl)
fj =fkj *flj
fjp =fkjp*flj
fjq =fkj *fik(2,ich(2),1,lj,y,min(2),bl)
fifj=fi*fj
dpp=0.0d0
dpq=0.0d0
dqp=0.0d0
dqq=0.0d0
if(dabs(G).gt.0.0d0)then
dpp=4.0d0*fip*fjp+qgpp*fifj+Gp*(fip*fj+fjp*fi)/G
dpq=4.0d0*fip*fjq+qgpq*fifj+(Gp*fjq*fi+Gq*fip*fj)/G
dqp=4.0d0*fiq*fjp+qgpq*fifj+(Gq*fjp*fi+Gp*fiq*fj)/G
dqq=4.0d0*fiq*fjq+qgqq*fifj+Gq*(fjq*fi+fiq*fj)/G
endif
t=0.125d0*(G11*dpp+G12*(dpq+dqp)+G22*dqq)*ds
v=en*fifj*ds
if(i.eq.2.and.j.eq.2)then
t22=t22+t
v22=v22+v
endif
if(i.eq.1.and.j.eq.1)then
t11=t11+t
v11=v11+v
endif
if(lh)then
h(i,j)=h(i,j)+t+v-emax2*fifj*ds
else
s(i,j)=s(i,j)+fi*fj*ds
h(i,j)=h(i,j)+t+v
endif
endif
2181 continue
endif
218 continue
do 219 i=1,nhdim
do 219 j=i+1,nhdim
s(i,j)=s(j,i)
219 h(i,j)=h(j,i)
endif
close(81)
if(lprop)close(816)
c
write(6,*)' H done',nhdim
write(6,*)t11*htocm,v11*htoc
write(6,*)nhdim,nhdim,h(nhdim,nhdim)*htocm,s(nhdim,nhdim)
write(6,*)
if(lwr)then
write(6,*)'Hamiltonian'
do 36 j=1,nhdim
36 write(6,601)(h(i,j)*htocm,i=1,nhdim)
endif
write(6,*)
if(kinddia.eq.0)then
write(6,*)' Basis not orthogonalized'
nfb=nhdim
CALL TRED12(nfb,h,c,q,2,IERR)
IF(IERR.NE.0)call report('diagonalization error')
else
write(6,*)' Basis orthogonalized'
CALL TRED12(nhdim,s,vu,d,2,IERR)
IF(IERR.NE.0)call report('diagonalization error')
write(6,*)'S-eigenvalues:'
write(6,601)(d(i),i=1,nhdim)
tolb=0.5d0
nfb=0
do 28 k=1,nhdim
if(abs(d(k)).gt.tolb)nfb=k
dk=dsqrt(d(k))
do 28 i=1,nhdim
28 vu(i,k)=vu(i,k)*dk
write(6,7007)nfb,tolb
7007 format(' Final dimension of the basis: ',i5,', tol: ',f12.6)
write(6,*)' Transforming H'
c
c First index trafo -use S as temporary matrix
do 29 k=1,nhdim
do 29 i=1,nhdim
ski=0.0d0
do 30 j=1,nhdim
30 ski=ski+vu(j,k)*h(j,i)
29 s(k,i)=ski
c
c Second index trafo
do 35 k=1,nhdim
do 35 l=1,nhdim
hkl=0.0d0
do 32 j=1,nhdim
32 hkl=hkl+vu(j,l)*s(k,j)
35 h(k,l)=hkl
write(6,*)' Diagonalizing'
CALL TRED12(nfb,h,dn,q,2,IERR)
IF(IERR.NE.0)call report('diagonalization error')
do 33 n=1,nfb
do 33 k=1,nhdim
cnk=0.0d0
do 34 i=1,nfb
34 cnk=cnk+dn(i,n)*vu(k,i)
33 c(k,n)=cnk
endif
if(lwr)then
write(6,*)
write(6,*)'Coefficients'
do 37 j=1,nhdim
37 write(6,601)(c(i,j),i=1,nhdim)
write(6,*)
endif
c
c calculate populations
qsum=0.0d0
akt=kelvin*3.16684d-6
do 703 i=1,nfb
703 qsum=qsum+exp(-q(i)/akt)
do 704 i=1,nfb
704 popul(i)=exp(-q(i)/akt)/qsum
open(8,file='ENERGIES.TXT')
do 702 i=1,nfb
702 write(8,723)i,q(i)*htocm,popul(i)
723 format(i4,f12.2,f12.5)
close(7)
write(6,*)'Energies (cm-1):'
write(6,601)(q(i)*htocm,i=1,nfb)
601 format(6f12.6)
write(6,*)'Transitions (cm-1):'
write(6,602)((q(i)-q(nfb))*htocm,i=1,nfb)
602 format(6f12.2)
c
c Wave function plots:
do 701 ii=1,ipl
ist=ip(ii)
write(6,*)ii,ist
write(s5,990)ist
990 format(i5)
do 65 is=1,5
65 if(s5(is:is).ne.' ')goto 651
651 open(81,file=s5(is:5)//'.PRN')
open(82,file=s5(is:5)//'b.PRN')
x=min(1)-dxh
xav=0.0d0
x2av=0.0d0
if(nvar.eq.1)then
do 8 ix=0,ngrh(1)
x=x+dxh
f=0.0d0
fb=fik(1,ich(1),0,ist,x,min(1),al)
do 9 k=1,nhdim
if(x.lt.0.0d0.and.x+dxh.ge.0.0d0)then
if(ist.eq.1)then
write(6,*)k,c(k,nfb-ist+1),fik(1,ich(1),0,k,x,min(1),al)
endif
endif
9 f=f+c(k,nfb-ist+1)*fik(1,ich(1),0,k,x,min(1),al)
xav =xav +x *f**2
x2av=x2av+x**2*f**2
write(81,8181)x,f
8 write(82,8181)x,fb
8181 format(2f12.6)
xav=xav*dxh
x2av=x2av*dxh
x2av=dsqrt(x2av-xav**2)
write(6,609)xav,x2av
609 format(' Average x: ',f12.5,' average dx: ',f12.5)
endif
if(nvar.eq.2)then
asum=0.0d0
yav=0.0d0
y2av=0.0d0
do 81 ix=0,ngrh(1)
x=x+dxh
if(ix.eq.0.or.ix.eq.ngrh(1))then
ddx=dxh/2.0d0
else
ddx=dxh
endif
y=min(2)-dyh
do 81 iy=0,ngrh(2)
y=y+dyh
if(iy.eq.0.or.iy.eq.ngrh(2))then
ddy=dyh/2.0d0
else
ddy=dyh
endif
ds=ddy*ddx
f=0.0d0
ih=0
do 91 k=1,nh(1)
fk=fik(1,ich(1),0,k,x,min(1),al)
do 91 l=1,nh(2)
ih=ih+1
91 f=f+c(ih,nfb-ist+1)*fk*fik(2,ich(2),0,l,y,min(2),bl)
ps=f**2*ds
if(ps.lt.1.0d-30)ps=0.0d0
xav =xav +x* ps
yav =yav +y* ps
y2av=y2av+y**2*ps
x2av=x2av+x**2*ps
asum=asum+ps
81 write(81,8182)x,y,f
8182 format(3G15.6)
write(6,*)'Sum: ',asum
x2av=dsqrt(x2av-xav**2)
y2av=dsqrt(y2av-yav**2)
write(6,6091)xav,yav,x2av,y2av
6091 format(' Average x,y: ',2f12.5,' average dx,dy: ',2f12.5)
endif
close(81)
close(82)
701 continue
c
c total probability plot:
xav=0.0d0
x2av=0.0d0
yav=0.0d0
y2av=0.0d0
it=int(kelvin)
write(s5,990)it
do 66 is=1,5
66 if(s5(is:is).ne.' ')goto 661
661 open(81,file=s5(is:5)//'P.PRN')
x=min(1)-dxh
plim=popul(nfb)*1.0d-4
if(nvar.eq.1)then
do 82 ix=0,ngrh(1)
x=x+dxh
f=0.0d0
do 92 kn=1,nfb
if(popul(kn).gt.plim)then
fifu=0.0d0
do 923 k=1,nhdim
923 fifu=fifu+c(k,kn)*fik(1,ich(1),0,k,x,min(1),al)
f=f+popul(kn)*fifu**2
endif
92 continue
xav=xav+x*f**2
x2av=x2av+x**2*f**2
82 write(81,8181)x,f
xav=xav*dxh
x2av=x2av*dxh
x2av=dsqrt(x2av-xav**2)
write(6,609)xav,x2av
endif
if(nvar.eq.2)then
do 812 ix=0,ngrh(1)
x=x+dxh
y=min(2)-dyh
do 812 iy=0,ngrh(2)
y=y+dyh
f=0.0d0
do 912 kn=1,nfb
if(popul(kn).gt.plim)then
fifu=0.0d0
ih=0
do 9124 k=1,nh(1)
fk=fik(1,ich(1),0,k,x,min(1),al)
do 9124 l=1,nh(2)
ih=ih+1
9124 fifu=fifu+c(ih,kn)*fk*fik(2,ich(2),0,l,y,min(2),bl)
f=f+popul(kn)*fifu**2
endif
912 continue
xav =xav +x* f*dxh*dyh
yav =yav +y* f*dxh*dyh
y2av=y2av+y**2*f*dxh*dyh
x2av=x2av+x**2*f*dxh*dyh
812 write(81,8182)x,y,f
x2av=dsqrt(x2av-xav**2)
y2av=dsqrt(y2av-yav**2)
write(6,6091)xav,yav,x2av,y2av
endif
close(81)
if(lprop)then
c
c Calculate for all states
write(6,*)' Calculating property'
do 51 ist=1,nfb
pa(ist)=0.0d0
51 pa2(ist)=0.0d0
open(816,file='PROPERTYNTH.LST')
x=min(1)-dxh
if(nvar.eq.1)then
do 52 ix=0,ngrh(1)
read(816,*)xdumm,prop
x=x+dxh
do 52 ist=1,nfb
f=0.0d0
do 53 k=1,nhdim
53 f=f+c(k,ist)*fik(1,ich(1),0,k,x,min(1),al)
pa(ist)=pa(ist)+f**2*prop*dxh
52 pa2(ist)=pa2(ist)+f**2*prop**2*dxh
endif
if(nvar.eq.2)then
do 54 ix=0,ngrh(1)
x=x+dxh
y=min(2)-dyh
do 54 iy=0,ngrh(2)
y=y+dyh
read(816,*)xdumm,ydumm,prop
do 54 ist=1,nfb
f=0.0d0
ih=0
do 55 k=1,nh(1)
fk=fik(1,ich(1),0,k,x,min(1),al)
do 55 l=1,nh(2)
ih=ih+1
55 f=f+c(ih,ist)*fk*fik(2,ich(2),0,l,y,min(2),bl)
pa(ist)=pa(ist)+f**2*prop*dxh*dyh
54 pa2(ist)=pa2(ist)+f**2*prop**2*dxh*dyh
endif
write(6,*)' :'
write(6,601)(pa(i),i=1,nfb)
write(6,*)' :'
write(6,601)(dsqrt(pa2(i)-pa(i)**2),i=1,nfb)
aave=0.0d0
dave=0.0d0
do 56 ist=1,nfb
aave=aave+popul(ist)*pa(ist)
56 dave=dave+popul(ist)*pa2(ist)
write(6,603)kelvin,aave,dsqrt(dave-aave**2)
603 format(' T = ',f7.1,', (T) = ',f12.6,', = ',f12.6)
endif
c
c time development of the wavefunction
if(nvar.eq.1)nh(2)=1
if(ltime)then
write(6,*)'Time propagation'
c
c initial packet
do 57 ist=1,nfb
f=0.0d0
if(nvar.eq.1)then
do 60 k=1,nh(1)
60 f=f+c(k,ist)*fik(1,ich(1),0,k,xt1,min(1),al)
endif
if(nvar.eq.2)then
ih=0
do 58 k=1,nh(1)
fk=fik(1,ich(1),0,k,xt1,min(1),al)
do 58 l=1,nh(2)
ih=ih+1
58 f=f+c(ih,ist)*fk*fik(2,ich(2),0,l,xt2,min(2),bl)
endif
dre(ist)=f
if(kelvin.gt.0.1d0)dre(ist)=dre(ist)*dsqrt(popul(ist))
57 continue
c
anorm=0.0d0
do 61 ist=1,nfb
61 anorm=anorm+dre(ist)**2
anorm=1.0d0/dsqrt(anorm)
do 62 ist=1,nfb
62 dre(ist)=dre(ist)*anorm
c
write(6,*)'Initial packet constructed'
c
t=0.0d0-dt
do 59 it=1,ntstep
t=t+dt
write(6,8009)it,t
8009 format(i4,' t = ',f22.4)
write(s5,990)it
do 64 is=1,5
64 if(s5(is:is).ne.' ')goto 641
641 open(81,file=s5(is:5)//'.t.prn')
if(nvar.eq.1)then
x=min(1)-dxh
do 63 ix=0,ngrh(1)
x=x+dxh
fr=0.0d0
fi=0.0d0
do 9123 ist=1,nfb
f=0.0d0
do 922 k=1,nhdim
922 f=f+c(k,nfb-ist+1)*fik(1,ich(1),0,k,x,min(1),al)
fr=fr+dre(ist)*f*cos(q(ist)*t)
9123 fi=fi+dre(ist)*f*sin(q(ist)*t)
63 write(81,8184)x,fr,fi,fr**2+fi**2
8184 format(4f12.6)
endif
if(nvar.eq.2)then
x=min(1)-dxh
do 811 ix=0,ngrh(1)
x=x+dxh
y=min(2)-dyh
do 811 iy=0,ngrh(2)
y=y+dyh
fr=0.0d0
fi=0.0d0
do 9121 ist=1,nfb
f=0.0d0
ih=0
do 911 k=1,nh(1)
fk=fik(1,ich(1),0,k,x,min(1),al)
do 911 l=1,nh(2)
ih=ih+1
911 f=f+c(ih,ist)*fk*fik(2,ich(2),0,l,y,min(2),bl)
fr=fr+dre(ist)*f*cos(q(ist)*t)
9121 fi=fi+dre(ist)*f*sin(q(ist)*t)
c wavefunction and probability at time t
811 write(81,8183)x,y,fr,fi,fr**2+fi**2
8183 format(5G15.6)
endif
59 close(81)
endif
stop
end
c ------------------------------------------------------------
subroutine report(s)
character*(*) s
write(6,*)s
stop
end
c ------------------------------------------------------------
function fik(ivar,ic,id,k,x,xs,al)
c
c basis function:
c ic ...0 plane waves
c 1 harmonic oscillator
c k: 1 2 3 4
c kindpw ... 0 cosinus const cos( pix/L) cos(2pix/L) cos(3pix/L)
c 1 cos and sin const sin(2pix/L) cos(2pix/L) sin(3pix/L)
c 2 sin sin(pix/L) sin(2pix/L) sin(3pix/L) sin(4pix/L)
c id=0 function
c id=1,2 first,second derivative
c k ...frequency (k=1 ... zero) or HO quantum number
c x ... variable
c al ... interval or HO frequency
c xs ... starting x for plane waves (origin)
c x0 ... cetral x for HO
implicit none
integer nh0,nvar0,kindpw,ivar
parameter (nh0=4100,nvar0=4)
real*8 fik,pi,x,al,sf,w,hk,hkp,hkpp,p,ank,y,xs,x0,am0,bohr,
1ff,dw
common/const/pi,bohr
common/ho/p(nh0,nh0),ank(nvar0,nh0),w(nvar0),x0(nvar0),am0(nvar0)
common/pw/kindpw
integer*4 ic,k,id,kk
if(ic.eq.0)then
c plane waves:
y=x-xs
if(k.eq.1.and.kindpw.ne.2)then
sf=1.0d0/dsqrt(al)
else
sf=dsqrt(2.0d0/al)
endif
if(id.eq.0)then
if(kindpw.eq.0)then
c K: 1 2 3 4
c const cos(pix/L) cos(2pix/L) cos(3pix/L) ...
fik=sf*cos(dble(k-1)*pi*y/al)
else
if(kindpw.eq.1)then
if(mod(k,2).ne.0)then
c K: 1 2 3 4 5 6
c const sin(2pix/L) cos(2pix/L) sin(4pix/L) cos(4pix/L) sin(6pix/L)
fik=sf*cos(dble(2*((k-1)/2))*pi*y/al)
else
fik=sf*sin(dble(2*(k/2))*pi*y/al)
endif
else
if(kindpw.eq.2)then
c K: 1 2 3 4
c sin(pix/L) sin(2pix/L) sin(3pix/L) sin(4pix/L)
fik=sf*sin(dble(k)*pi*y/al)
else
write(6,*)
call report('unknown function requested')
endif
endif
endif
else
if(id.eq.1)then
if(kindpw.eq.0)then
fik=-sf*sin(dble(k-1)*pi*y/al)*dble(k-1)/al*pi
else
if(kindpw.eq.1)then
if(mod(k,2).ne.0)then
fik=-sf*sin(dble(2*((k-1)/2))*pi*y/al)
1 *dble(2*((k-1)/2))*pi/al
else
fik= sf*cos(dble(2*(k/2))*pi*y/al)*dble(2*(k/2))*pi/al
endif
else
fik=sf*cos(dble(k)*pi*y/al)*dble(k)/al*pi
endif
endif
else
if(id.eq.2)then
if(kindpw.eq.0)then
fik=-sf*cos(dble(k-1)*pi*y/al)*(dble(k-1)*pi/al)**2
else
if(kindpw.eq.1)then
if(mod(k,2).ne.0)then
fik=-sf*cos(dble(2*((k-1)/2))*pi*y/al)
1 *(dble(2*((k-1)/2))*pi/al)**2
else
fik=-sf*sin(dble(2*(k/2))*pi*y/al)*(dble(2*(k/2))*pi/al)**2
endif
else
fik=-sf*sin(dble(k)*pi*y/al)*(dble(k)*pi/al)**2
endif
endif
else
call report('higher ders not implemented')
endif
endif
endif
else
dw=dsqrt(w(ivar)*am0(ivar))/bohr
y=dw*(x-x0(ivar))
ff=ank(ivar,k)*exp(-y**2/2.0d0)
hk=0.0d0
hkp=0.0d0
hkpp=0.0d0
do 1 kk=1,k
1 hk=hk+p(k,kk)*y**(kk-1)
if(id.gt.0)then
do 2 kk=2,k
2 hkp=hkp+dble(kk-1)*p(k,kk)*y**(kk-2)
if(id.gt.1)then
do 3 kk=3,k
3 hkpp=hkpp+dble((kk-1)*(kk-2))*p(k,kk)*y**(kk-3)
endif
endif
if(id.eq.0)then
fik=ff*hk
else
if(id.eq.1)then
fik=ff*(hkp-y*hk)*dw
else
if(id.eq.2)then
fik=ff*(hk*(y**2-1.0d0)-2.0d0*y*hkp+hkpp)*dw**2
else
call report('higher ders not implemented')
endif
endif
endif
endif
return
end
c ------------------------------------------------------------
SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR)
IMPLICIT REAL*8(A-H,O-Z)
IMPLICIT INTEGER*4 (I-N)
PARAMETER (nhdim0=4100)
DIMENSION A(nhdim0,N),Z(nhdim0,N),D(nhdim0),E(nhdim0)
COMMON/DIAG/E
IF (IEIG .LT. 0) GO TO 110
DO 100 I = 1, N
C
DO 100 J = 1, I
Z(I,J) = A(I,J)
100 CONTINUE
C
110 ICODE = IABS(IEIG)
IF (N .EQ. 1) GO TO 320
C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- ::::::::::
DO 300 II = 2, N
I = N + 2 - II
L = I - 1
H = 0.0D0
SCALE = 0.0D0
IF (L .LT. 2) GO TO 130
C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) ::::::::::
DO 120 K = 1, L
120 SCALE = SCALE + DABS(Z(I,K))
C
IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140
130 E(I) = Z(I,L)
GO TO 290
C
140 DO 150 K = 1, L
F = Z(I,K) / SCALE
Z(I,K) = F
H = H + F * F
150 CONTINUE
C
G = -DSIGN(DSQRT(H),F)
E(I) = SCALE * G
H = H - F * G
Z(I,L) = F - G
F = 0.0D0
C
DO 240 J = 1, L
IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H)
G = 0.0D0
C :::::::::: FORM ELEMENT OF A*U ::::::::::
DO 180 K = 1, J
180 G = G + Z(J,K) * Z(I,K)
C
JP1 = J + 1
IF (L .LT. JP1) GO TO 220
C
DO 200 K = JP1, L
200 G = G + Z(K,J) * Z(I,K)
C :::::::::: FORM ELEMENT OF P ::::::::::
220 E(J) = G / H
F = F + E(J) * Z(I,J)
240 CONTINUE
C
HH = F / (H + H)
C :::::::::: FORM REDUCED A ::::::::::
DO 260 J = 1, L
F = Z(I,J)
G = E(J) - HH * F
E(J) = G
C
DO 260 K = 1, J
Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
260 CONTINUE
C
DO 280 K = 1, L
280 Z(I,K) = SCALE * Z(I,K)
C
290 D(I) = H
300 CONTINUE
C
320 D(1) = 0.0D0
E(1) = 0.0D0
C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES ::::::::::
IF (ICODE .NE. 2) GO TO 600
DO 500 I = 1, N
L = I - 1
IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380
C
DO 360 J = 1, L
G = 0.0D0
C
DO 340 K = 1, L
340 G = G + Z(I,K) * Z(K,J)
C
DO 360 K = 1, L
Z(K,J) = Z(K,J) - G * Z(K,I)
360 CONTINUE
C
380 D(I) = Z(I,I)
Z(I,I) = 1.0D0
IF (L .LT. 1) GO TO 500
C
DO 400 J = 1, L
Z(I,J) = 0.0D0
Z(J,I) = 0.0D0
400 CONTINUE
C
500 CONTINUE
C
620 CALL TQL12 (N,Z,D,IERR,ICODE)
RETURN
C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY ::::::::::
600 DO 610 I=1,N
610 D(I) = Z(I,I)
GO TO 620
END
C
SUBROUTINE TQL12(N,Z,D,IERR,ICODE)
IMPLICIT REAL*8(A-H,O-Z)
IMPLICIT INTEGER*4 (I-N)
PARAMETER (nhdim0=4100)
DIMENSION Z(nhdim0,N),D(nhdim0),E(nhdim0)
REAL*8 MACHEP
COMMON/DIAG/E
EPS = 1.0D0
10 EPS = 0.50D0*EPS
TOL1 = EPS + 1.0D0
IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10
IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS
MACHEP = EPS
C
C MACHEP=16.0D0**(-13)
C
IERR = 0
IF (N .EQ. 1) GO TO 1001
C
DO 100 I = 2, N
100 E(I-1) = E(I)
C
F = 0.0D0
B = 0.0D0
E(N) = 0.0D0
C
DO 240 L = 1, N
J = 0
H = MACHEP * (DABS(D(L)) + DABS(E(L)))
IF (B .LT. H) B = H
C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT ::::::::::
DO 110 M = L, N
IF (DABS(E(M)) .LE. B) GO TO 120
C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C THROUGH THE BOTTOM OF THE LOOP ::::::::::
110 CONTINUE
C
120 IF (M .EQ. L) GO TO 220
130 IF (J .EQ. 30) GO TO 1000
J = J + 1
C :::::::::: FORM SHIFT ::::::::::
L1 = L + 1
G = D(L)
P = (D(L1) - G) / (2.0D0 * E(L))
R = DSQRT(P*P+1.0D0)
D(L) = E(L) / (P + DSIGN(R,P))
H = G - D(L)
C
DO 140 I = L1, N
140 D(I) = D(I) - H
C
F = F + H
C :::::::::: QL TRANSFORMATION ::::::::::
P = D(M)
C = 1.0D0
S = 0.0D0
MML = M - L
C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- ::::::::::
DO 200 II = 1, MML
I = M - II
G = C * E(I)
H = C * P
IF (DABS(P) .LT. DABS(E(I))) GO TO 150
C = E(I) / P
R = DSQRT(C*C+1.0D0)
E(I+1) = S * P * R
S = C / R
C = 1.0D0 / R
GO TO 160
150 C = P / E(I)
R = DSQRT(C*C+1.0D0)
E(I+1) = S * E(I) * R
S = 1.0D0 / R
C = C * S
160 P = C * D(I) - S * G
D(I+1) = H + S * (C * G + S * D(I))
C :::::::::: FORM VECTOR ::::::::::
IF (ICODE .NE. 2) GO TO 200
DO 180 K = 1, N
H = Z(K,I+1)
Z(K,I+1) = S * Z(K,I) + C * H
Z(K,I) = C * Z(K,I) - S * H
180 CONTINUE
C
200 CONTINUE
C
E(L) = S * P
D(L) = C * P
IF (DABS(E(L)) .GT. B) GO TO 130
220 D(L) = D(L) + F
240 CONTINUE
C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS ::::::::::
DO 300 II = 2, N
I = II - 1
K = I
P = D(I)
C
DO 260 J = II, N
IF (D(J) .LE. P) GO TO 260
K = J
P = D(J)
260 CONTINUE
C
IF (K .EQ. I) GO TO 300
D(K) = D(I)
D(I) = P
C
IF (ICODE .NE. 2) GO TO 300
DO 280 J = 1, N
P = Z(J,I)
Z(J,I) = Z(J,K)
Z(J,K) = P
280 CONTINUE
C
300 CONTINUE
C
GO TO 1001
C :::::::::: SET ERROR -- NO CONVERGENCE TO AN
C EIGENVALUE AFTER 30 ITERATIONS ::::::::::
1000 IERR = L
1001 RETURN
END
c ------------------------------------------------------------------
subroutine order(d1,d2,d3,i1,i2,i3)
real*8 d1,d2,d3,t
integer*4 i1,i2,i3
if(d1.gt.d2)then
t=d1
d1=d2
d2=t
it=i1
i1=i2
i2=it
endif
if(d1.gt.d3)then
t=d1
d1=d3
d3=t
it=i1
i1=i3
i3=it
endif
if(d2.gt.d3)then
t=d2
d2=d3
d3=t
it=i2
i2=i3
i3=it
endif
return
end
c ------------------------------------------------------------
subroutine abc(a,b,c,q1,q2,q3,x1,x2,x3)
real*8 a,b,c,q1,q2,q3,x1,x2,x3,down
down=-(q1-q2)*(q2-q3)*(q3-q1)
a=(x1*(q2-q3)+x2*(q3-q1)+x3*(q1-q2))/down
b=-(x1*(q2**2-q3**2)+x2*(q3**2-q1**2)+x3*(q1**2-q2**2))/down
c=(x1*(q2-q3)*q2*q3+x2*(q3-q1)*q1*q3+x3*(q1-q2)*q1*q2)/down
return
end
c ------------------------------------------------------------
subroutine find(i1,i2,i3,var,x,lp1,ie,tp)
c finds 3 y-arrays in x1 x2 x3 closest to x
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
logical lp1,nonequal
real*8 var(*)
to=1.0d-4
c
c initialize x1 as var(1)
i1=1
d1=ddi(lp1,var(i1),x,tp)
c
c get some x2 different from x1
do 10 i2=2,ie
10 if(nonequal(lp1,var(i2),var(i1),to,tp))goto 11
call report('x2 not found in findlist')
11 d2=ddi(lp1,var(i2),x,tp)
c
c get x3 different from x1 and x2
do 12 i3=3,ie
12 if(nonequal(lp1,var(i3),var(i1),to,tp).and.
1nonequal(lp1,var(i3),var(i2),to,tp))goto 13
call report('x3 not found in findlist')
13 d3=ddi(lp1,var(i3),x,tp)
call order(d1,d2,d3,i1,i2,i3)
c look if some points are closer:
do 1 ii=4,ie
dist=ddi(lp1,var(ii),x,tp)
if(dist.lt.d1-to)then
d3=d2
i3=i2
d2=d1
i2=i1
d1=dist
i1=ii
endif
if(dist.lt.d2-to.and.nonequal(lp1,var(ii),var(i1),to,tp))then
d3=d2
i3=i2
d2=dist
i2=ii
endif
if(dist.lt.d3-to.and.nonequal(lp1,var(ii),var(i1),to,tp).and.
1nonequal(lp1,var(ii),var(i2),to,tp))then
d3=dist
i3=ii
endif
1 continue
return
end
c ------------------------------------------------------------
subroutine findx(i1,i2,i3,var,x,lp1,ie,ngp0,nvar0,tp)
c finds 3 y-arrays in x1 x2 x3 closest to x
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
logical lp1,nonequal
real*8 var(ngp0,nvar0)
to=1.0d-4
c
c initialize x1 as var(1
i1=1
d1=ddi(lp1,var(i1,1),x,tp)
c
c get some x2 different from x1
do 10 i2=2,ie
10 if(nonequal(lp1,var(i2,1),var(i1,1),to,tp))goto 11
call report('x2 not found in findx')
11 d2=ddi(lp1,var(i2,1),x,tp)
c
c get x3 different from x1 and x2
do 12 i3=3,ie
12 if(nonequal(lp1,var(i3,1),var(i1,1),to,tp).and.
1nonequal(lp1,var(i3,1),var(i2,1),to,tp))goto 13
call report('x3 not found in findx')
13 d3=ddi(lp1,var(i3,1),x,tp)
call order(d1,d2,d3,i1,i2,i3)
c look if some points are closer:
do 1 ii=4,ie
dist=ddi(lp1,var(ii,1),x,tp)
if(dist.lt.d1-to)then
d3=d2
i3=i2
d2=d1
i2=i1
d1=dist
i1=ii
endif
if(dist.lt.d2-to.and.nonequal(lp1,var(ii,1),var(i1,1),to,tp))then
d3=d2
i3=i2
d2=dist
i2=ii
endif
if(dist.lt.d3-to.and.nonequal(lp1,var(ii,1),var(i1,1),to,tp).and.
1nonequal(lp1,var(ii,1),var(i2,1),to,tp))then
d3=dist
i3=ii
endif
1 continue
return
end
c ------------------------------------------------------------
subroutine findlist(l1,l2,l3,var,x,lp1,ie,ngp0,nvar0,tp)
c finds 3 y-arrays in x1 x2 x3 closest to x
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
logical lp1,nonequal
real*8 var(ngp0,nvar0)
dimension l1(*),l2(*),l3(*)
to=1.0d-4
call findx(i1,i2,i3,var,x,lp1,ie,ngp0,nvar0,tp)
c record list of ys for xs same as x1,x2,x3
nl1=1
nl2=1
nl3=1
do 2 ii=1,ie
if(.not.nonequal(lp1,var(ii,1),var(i1,1),to,tp))then
nl1=nl1+1
l1(nl1)=ii
endif
if(.not.nonequal(lp1,var(ii,1),var(i2,1),to,tp))then
nl2=nl2+1
l2(nl2)=ii
endif
if(.not.nonequal(lp1,var(ii,1),var(i3,1),to,tp))then
nl3=nl3+1
l3(nl3)=ii
endif
2 continue
l1(1)=nl1
l2(1)=nl2
l3(1)=nl3
return
end
c ------------------------------------------------------------
subroutine findy(y,l1,i1,i2,i3,var,lp,ngp0,nvar0,tp)
c in l1 find y1 y2 y3 closest to y
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
logical lp,nonequal
real*8 var(ngp0,nvar0)
dimension l1(*)
to=1.0d-4
c initalize y1
i1=2
d1=ddi(lp,var(l1(i1),2),y,tp)
c find y2 different from y1
do 10 i2=3,l1(1)
10 if(nonequal(lp,var(l1(i2),2),var(l1(i1),2),to,tp))goto 11
call report('y2 not found in findy')
11 d2=ddi(lp,var(l1(i2),2),y,tp)
c find y3 different from y1 and from y2
do 12 i3=4,l1(1)
12 if(nonequal(lp,var(l1(i3),2),var(l1(i1),2),to,tp).and.
1nonequal(lp,var(l1(i3),2),var(l1(i2),2),to,tp))goto 13
call report('y3 not found in findy')
13 d3=ddi(lp,var(l1(i3),2),y,tp)
call order(d1,d2,d3,i1,i2,i3)
c look for closer points:
do 1 ii=5,l1(1)
dist=ddi(lp,var(l1(ii),2),y,tp)
if(dist.lt.d1-to)then
d3=d2
i3=i2
d2=d1
i2=i1
d1=dist
i1=ii
endif
if(dist.lt.d2-to.and.nonequal(lp,var(l1(ii),2),var(l1(i1),2),to,tp
1))then
d3=d2
i3=i2
d2=dist
i2=ii
endif
if(dist.lt.d3-to.and.nonequal(lp,var(l1(ii),2),var(l1(i1),2),to,tp
1).and.nonequal(lp,var(l1(ii),2),var(l1(i2),2),to,tp))then
d3=dist
i3=ii
endif
1 continue
return
end
c ---------------------------------------------------------------
function ddi(l,x,y,tp)
real*8 x,y,tp,d,ddi
logical l
d=abs(x-y)
if(l)then
if(d.gt.abs(x-y-tp))d=abs(x-y-tp)
if(d.gt.abs(x-y+tp))d=abs(x-y+tp)
endif
ddi=d
return
end
c ------------------------------------------------------------
function nonequal(l,x,y,to,tp)
logical nonequal,l
real*8 x,y,to,tp
if((l.and.abs(x-y).gt.to.and.abs(x-y+tp).gt.to.and.
1abs(x-y-tp).gt.to).or.(.not.l.and.abs(x-y).gt.to))then
nonequal=.true.
else
nonequal=.false.
endif
return
end
c ------------------------------------------------------------
subroutine getpar(nvar,nvar0,lper,ngr,min,max,nbase,nbas,nh,ngrh,
1ice,ic,ich,wcm,x0,am0,ipl,ip,ikin,kindpw,kinddia,T,w,nbas0,nh0,
2lh,emax1,emax2,peri,ltime,lwr)
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
character*4 key
c icon=0 ... reads par for surfaces
c icon=1 ... reads par for peswf
dimension ip(*),ngr(*),min(*),max(*),ngrh(*),nbase(*),
1nbas(*),nh(*),ice(*),ic(*),ich(*),x0(*),am0(*),w(*),lper(*),
2peri(*)
logical lper,lh,ltime,lwr
real*8 min,max
common/timec/dt,xt1,xt2,ntstep
lwr=.false.
ipl=0
ntstep=10
xt1=0.0d0
xt2=0.0d0
dt=4.30d0
lh=.false.
ikin=1
kindpw=1
kinddia=1
T=298.15d0
ltime=.false.
il=0
open(8,file='SURF.OPT')
1 read(8,800,end=88,err=88)key
800 format(a4)
il=il+1
c
c Obligatory Parameters:
if(key.eq.'NVAR'.and.il.eq.1)then
read(8,*)nvar
if(nvar.gt.nvar0)call report('Too many variables')
do 2 ii=1,nvar
write(6,6005)ii
6005 format(' Variable',i2,':')
peri(ii)=360.0d0
read(8,*)lper(ii)
read(8,*)ngr(ii)
read(8,*)min(ii),max(ii)
read(8,*)ngrh(ii)
read(8,*)nbase(ii),nbas(ii),nh(ii)
read(8,*)ice(ii),ic(ii),ich(ii)
read(8,*)wcm,x0(ii),am0(ii)
if(ice(ii).eq.0)then
write(6,6002)nbase(ii), ' plane waves','E'
6002 format(i4,a12,' for ',a1)
else
if(ice(ii).eq.1)then
write(6,6002)nbase(ii),' harmonics','E'
else
write(6,*)'Quadratic fit for E'
endif
endif
if(ic(ii).eq.0)then
write(6,6002)nbas(ii), ' plane waves','G'
else
if(ic(ii).eq.1)then
write(6,6002)nbas(ii),' harmonics','G'
else
write(6,*)'Quadratic fit for G'
endif
endif
if(ich(ii).eq.0)then
write(6,6002)nh(ii), ' plane waves','H'
else
if(ich(ii).eq.1)then
write(6,6002)nh(ii),' harmonics','H'
else
call report('Unknown basis option for H')
endif
endif
if(nbase(ii).gt.nbas0)call report('too many E basis functions')
if(nbas(ii). gt.nbas0)call report('too many G basis functions')
if(nh(ii). gt.nh0 )call report('too many H basis functions')
if(ic(ii).eq.1.or.ich(ii).eq.1)write(6,5006)wcm,x0(ii),am0(ii)
5006 format(' HO basis',f12.2,' cm-1, x0 = ',f12.6,' mass = ',f10.2)
2 w(ii)=wcm/htocm
else
if(il.eq.1)call report('NVAR must be defined at the first line.')
endif
c
c Optional Parameters
if(key.eq.'STAT')read(8,*)ipl,(ip(ii),ii=1,ipl)
if(key.eq.'KINE')read(8,*)ikin
if(key.eq.'PLAN')read(8,*)kindpw
if(key.eq.'DIAG')read(8,*)kinddia
if(key.eq.'TEMP')read(8,*)T
if(key.eq.'TDEP')read(8,*)ltime
if(key.eq.'XT10')read(8,*)xt1
if(key.eq.'XT20')read(8,*)xt2
if(key.eq.'NTST')read(8,*)ntstep
if(key.eq.'DTST')read(8,*)dt
if(key.eq.'PERI')read(8,*)(peri(i),i=1,NVAR)
c large output:
if(key.eq.'LWRI')read(8,*)lwr
if(key.eq.'HOLE')then
lh=.true.
read(8,*)emax1,emax2
endif
goto 1
88 close(8)
return
end
c ------------------------------------------------------------
subroutine check(ipl,ikin,nvar,ic,kindpw)
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
dimension ic(*)
write(6,*)ipl,' states will be plotted'
if(ikin.eq.2)then
write(6,*)' Normal kinetic operator'
else
if(ikin.eq.1)then
write(6,*)' Kinetic operator per partes'
else
call report('Unknown T')
endif
endif
if(ikin.eq.2.and.nvar.gt.1)
1call report('Multidimensions and normal kinetic not implemented')
if(nvar.gt.1.and.(ic(1).ne.0.and.ic(1).ne.1))
1call report('With multidimensions fit required for G')
if(kindpw.eq.0)then
write(6,*)' Cosinus PW'
else
if(kindpw.eq.1)then
write(6,*)' Cosinus and sinus PW'
else
if(kindpw.eq.2)then
write(6,*)' Sinus PW'
else
call report('Unknown PW option')
endif
endif
endif
return
end