9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00

Removed labels in loc.f

This commit is contained in:
Anthony Scemama 2022-11-19 16:16:04 +01:00
parent 8245c7b3f7
commit 2798183c6d
2 changed files with 163 additions and 142 deletions

View File

@ -1095,9 +1095,9 @@ double precision function overlap_orb_ylm_grid(nptsgrid,r_orb,npower_orb,center_
implicit none implicit none
!! PSEUDOS !! PSEUDOS
integer nptsgridmax,nptsgrid integer nptsgridmax,nptsgrid
double precision coefs_pseudo,ptsgrid
parameter(nptsgridmax=50) parameter(nptsgridmax=50)
common/pseudos/coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3) double precision coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3)
common/pseudos/coefs_pseudo,ptsgrid
!!!!! !!!!!
integer npower_orb(3),l,m,i integer npower_orb(3),l,m,i
double precision x,g_orb,two_pi,dx,dphi,term,orb_phi,ylm_real,sintheta,r_orb,phi,center_orb(3) double precision x,g_orb,two_pi,dx,dphi,term,orb_phi,ylm_real,sintheta,r_orb,phi,center_orb(3)
@ -1235,10 +1235,10 @@ end
subroutine initpseudos(nptsgrid) subroutine initpseudos(nptsgrid)
implicit none implicit none
integer nptsgridmax,nptsgrid,ik integer nptsgridmax,nptsgrid,ik
double precision coefs_pseudo,ptsgrid
double precision p,q,r,s double precision p,q,r,s
parameter(nptsgridmax=50) parameter(nptsgridmax=50)
common/pseudos/coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3) double precision :: coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3)
common/pseudos/coefs_pseudo,ptsgrid
p=1.d0/dsqrt(2.d0) p=1.d0/dsqrt(2.d0)
q=1.d0/dsqrt(3.d0) q=1.d0/dsqrt(3.d0)

View File

@ -1,6 +1,6 @@
c************************************************************************ c************************************************************************
subroutine maxovl(n,m,s,t,w) subroutine maxovl(n,m,s,t,w)
C C
C This subprogram contains an iterative procedure to find the C This subprogram contains an iterative procedure to find the
C unitary transformation of a set of n vectors which maximizes C unitary transformation of a set of n vectors which maximizes
C the sum of their square overlaps with a set of m reference C the sum of their square overlaps with a set of m reference
@ -10,7 +10,7 @@ C S: overlap matrix <ref|vec>
C T: rotation matrix C T: rotation matrix
C W: new overlap matrix C W: new overlap matrix
C C
C C
implicit real*8(a-h,o-y),logical*1(z) implicit real*8(a-h,o-y),logical*1(z)
! parameter (id1=700) ! parameter (id1=700)
! dimension s(id1,id1),t(id1,id1),w(id1,id1) ! dimension s(id1,id1),t(id1,id1),w(id1,id1)
@ -29,23 +29,26 @@ C conv=1.d-6
* 5x,'following the principle of maximum overlap with a set of', * 5x,'following the principle of maximum overlap with a set of',
* i3,' reference vectors'/5x,'required convergence on rotation ', * i3,' reference vectors'/5x,'required convergence on rotation ',
* 'angle =',f13.10///5x,'Starting overlap matrix'/) * 'angle =',f13.10///5x,'Starting overlap matrix'/)
do 6 i=1,m do i=1,m
write (6,145) i write (6,145) i
6 write (6,150) (s(i,j),j=1,n) write (6,150) (s(i,j),j=1,n)
end do
8 mm=m-1 8 mm=m-1
if (m.lt.n) mm=m if (m.lt.n) mm=m
iter=0 iter=0
do 20 j=1,n do j=1,n
do 16 i=1,n do i=1,n
t(i,j)=0.d0 t(i,j)=0.d0
16 continue end do
do 18 i=1,m do i=1,m
18 w(i,j)=s(i,j) w(i,j)=s(i,j)
20 t(j,j)=1.d0 enddo
t(j,j)=1.d0
enddo
sum=0.d0 sum=0.d0
do 10 i=1,m do i=1,m
sum=sum+s(i,i)*s(i,i) sum=sum+s(i,i)*s(i,i)
10 continue end do
sum=sum/m sum=sum/m
if (zprt) write (6,12) sum if (zprt) write (6,12) sum
12 format (//5x,'Average square overlap =',f10.6) 12 format (//5x,'Average square overlap =',f10.6)
@ -54,18 +57,18 @@ C conv=1.d-6
j=1 j=1
21 if (j.ge.last) goto 30 21 if (j.ge.last) goto 30
sum=0.d0 sum=0.d0
do i=1,n
do 22 i=1,n sum=sum+s(i,j)*s(i,j)
22 sum=sum+s(i,j)*s(i,j) enddo
if (sum.gt.small) goto 28 if (sum.gt.small) goto 28
do 24 i=1,n do i=1,n
sij=s(i,j) sij=s(i,j)
s(i,j)=-s(i,last) s(i,j)=-s(i,last)
s(i,last)=sij s(i,last)=sij
tij=t(i,j) tij=t(i,j)
t(i,j)=-t(i,last) t(i,j)=-t(i,last)
t(i,last)=tij t(i,last)=tij
24 continue end do
last=last-1 last=last-1
goto 21 goto 21
28 j=j+1 28 j=j+1
@ -101,17 +104,18 @@ C conv=1.d-6
sine=1.d0 sine=1.d0
34 delta=sine*(a*sine+b*cosine) 34 delta=sine*(a*sine+b*cosine)
if (zprt.and.delta.lt.0.d0) write (6,71) i,j,a,b,sine,cosine,delta if (zprt.and.delta.lt.0.d0) write (6,71) i,j,a,b,sine,cosine,delta
do 35 k=1,m do k=1,m
p=s(k,i)*cosine-s(k,j)*sine p=s(k,i)*cosine-s(k,j)*sine
q=s(k,i)*sine+s(k,j)*cosine q=s(k,i)*sine+s(k,j)*cosine
s(k,i)=p s(k,i)=p
35 s(k,j)=q s(k,j)=q
do 40 k=1,n enddo
p=t(k,i)*cosine-t(k,j)*sine do k=1,n
q=t(k,i)*sine+t(k,j)*cosine p=t(k,i)*cosine-t(k,j)*sine
t(k,i)=p q=t(k,i)*sine+t(k,j)*cosine
t(k,j)=q t(k,i)=p
40 continue t(k,j)=q
enddo
45 d=dabs(sine) 45 d=dabs(sine)
if (d.le.amax) goto 50 if (d.le.amax) goto 50
imax=i imax=i
@ -132,43 +136,50 @@ C conv=1.d-6
* 'in subroutine maxovl ***'//) * 'in subroutine maxovl ***'//)
stop stop
100 continue 100 continue
do 120 j=1,n do j=1,n
if (s(j,j).gt.0.d0) goto 120 if (s(j,j).gt.0.d0) cycle
do 105 i=1,m do i=1,m
105 s(i,j)=-s(i,j) s(i,j)=-s(i,j)
do 110 i=1,n enddo
110 t(i,j)=-t(i,j) do i=1,n
120 continue t(i,j)=-t(i,j)
enddo
enddo
sum=0.d0 sum=0.d0
do 125 i=1,m do i=1,m
125 sum=sum+s(i,i)*s(i,i) sum=sum+s(i,i)*s(i,i)
enddo
sum=sum/m sum=sum/m
do 122 i=1,m do i=1,m
do 122 j=1,n do j=1,n
sw=s(i,j) sw=s(i,j)
s(i,j)=w(i,j) s(i,j)=w(i,j)
122 w(i,j)=sw w(i,j)=sw
enddo
enddo
if (.not.zprt) return if (.not.zprt) return
write (6,12) sum write (6,12) sum
write (6,130) write (6,130)
130 format (//5x,'transformation matrix') 130 format (//5x,'transformation matrix')
do 140 i=1,n do i=1,n
write (6,145) i write (6,145) i
140 write (6,150) (t(i,j),j=1,n) write (6,150) (t(i,j),j=1,n)
enddo
145 format (i8) 145 format (i8)
150 format (2x,10f12.8) 150 format (2x,10f12.8)
write (6,160) write (6,160)
160 format (//5x,'new overlap matrix'/) 160 format (//5x,'new overlap matrix'/)
do 170 i=1,m do i=1,m
write (6,145) i write (6,145) i
170 write (6,150) (w(i,j),j=1,n) write (6,150) (w(i,j),j=1,n)
enddo
return return
end end
c************************************************************************ c************************************************************************
subroutine maxovl_no_print(n,m,s,t,w) subroutine maxovl_no_print(n,m,s,t,w)
C C
C This subprogram contains an iterative procedure to find the C This subprogram contains an iterative procedure to find the
C unitary transformation of a set of n vectors which maximizes C unitary transformation of a set of n vectors which maximizes
C the sum of their square overlaps with a set of m reference C the sum of their square overlaps with a set of m reference
@ -178,7 +189,7 @@ C S: overlap matrix <ref|vec>
C T: rotation matrix C T: rotation matrix
C W: new overlap matrix C W: new overlap matrix
C C
C C
implicit real*8(a-h,o-y),logical*1(z) implicit real*8(a-h,o-y),logical*1(z)
parameter (id1=300) parameter (id1=300)
dimension s(id1,id1),t(id1,id1),w(id1,id1) dimension s(id1,id1),t(id1,id1),w(id1,id1)
@ -193,17 +204,19 @@ C conv=1.d-6
8 mm=m-1 8 mm=m-1
if (m.lt.n) mm=m if (m.lt.n) mm=m
iter=0 iter=0
do 20 j=1,n do j=1,n
do 16 i=1,n do i=1,n
t(i,j)=0.d0 t(i,j)=0.d0
16 continue enddo
do 18 i=1,m do i=1,m
18 w(i,j)=s(i,j) w(i,j)=s(i,j)
20 t(j,j)=1.d0 enddo
t(j,j)=1.d0
enddo
sum=0.d0 sum=0.d0
do 10 i=1,m do i=1,m
sum=sum+s(i,i)*s(i,i) sum=sum+s(i,i)*s(i,i)
10 continue enddo
sum=sum/m sum=sum/m
12 format (//5x,'Average square overlap =',f10.6) 12 format (//5x,'Average square overlap =',f10.6)
if (n.eq.1) goto 100 if (n.eq.1) goto 100
@ -211,18 +224,19 @@ C conv=1.d-6
j=1 j=1
21 if (j.ge.last) goto 30 21 if (j.ge.last) goto 30
sum=0.d0 sum=0.d0
do 22 i=1,n do i=1,n
22 sum=sum+s(i,j)*s(i,j) sum=sum+s(i,j)*s(i,j)
enddo
if (sum.gt.small) goto 28 if (sum.gt.small) goto 28
do 24 i=1,n do i=1,n
sij=s(i,j) sij=s(i,j)
s(i,j)=-s(i,last) s(i,j)=-s(i,last)
s(i,last)=sij s(i,last)=sij
tij=t(i,j) tij=t(i,j)
t(i,j)=-t(i,last) t(i,j)=-t(i,last)
t(i,last)=tij t(i,last)=tij
24 continue end do
last=last-1 last=last-1
goto 21 goto 21
28 j=j+1 28 j=j+1
@ -232,50 +246,52 @@ C conv=1.d-6
jmax=0 jmax=0
dmax=0.d0 dmax=0.d0
amax=0.d0 amax=0.d0
do 60 i=1,mm do i=1,mm
ip=i+1 ip=i+1
do 50 j=ip,n do j=ip,n
a=s(i,j)*s(i,j)-s(i,i)*s(i,i) a=s(i,j)*s(i,j)-s(i,i)*s(i,i)
b=-s(i,i)*s(i,j) b=-s(i,i)*s(i,j)
if (j.gt.m) goto 31 if (j.gt.m) goto 31
a=a+s(j,i)*s(j,i)-s(j,j)*s(j,j) a=a+s(j,i)*s(j,i)-s(j,j)*s(j,j)
b=b+s(j,i)*s(j,j) b=b+s(j,i)*s(j,j)
31 b=b+b 31 b=b+b
if (a.eq.0.d0) goto 32 if (a.eq.0.d0) goto 32
ba=b/a ba=b/a
if (dabs(ba).gt.small) goto 32 if (dabs(ba).gt.small) goto 32
if (a.gt.0.d0) goto 33 if (a.gt.0.d0) goto 33
tang=-0.5d0*ba tang=-0.5d0*ba
cosine=1.d0/dsqrt(1.d0+tang*tang) cosine=1.d0/dsqrt(1.d0+tang*tang)
sine=tang*cosine sine=tang*cosine
goto 34 goto 34
32 tang=0.d0 32 tang=0.d0
if (b.ne.0.d0) tang=(a+dsqrt(a*a+b*b))/b if (b.ne.0.d0) tang=(a+dsqrt(a*a+b*b))/b
cosine=1.d0/dsqrt(1.d0+tang*tang) cosine=1.d0/dsqrt(1.d0+tang*tang)
sine=tang*cosine sine=tang*cosine
goto 34 goto 34
33 cosine=0.d0 33 cosine=0.d0
sine=1.d0 sine=1.d0
34 delta=sine*(a*sine+b*cosine) 34 delta=sine*(a*sine+b*cosine)
do 35 k=1,m do k=1,m
p=s(k,i)*cosine-s(k,j)*sine p=s(k,i)*cosine-s(k,j)*sine
q=s(k,i)*sine+s(k,j)*cosine q=s(k,i)*sine+s(k,j)*cosine
s(k,i)=p s(k,i)=p
35 s(k,j)=q s(k,j)=q
do 40 k=1,n enddo
p=t(k,i)*cosine-t(k,j)*sine do k=1,n
q=t(k,i)*sine+t(k,j)*cosine p=t(k,i)*cosine-t(k,j)*sine
t(k,i)=p q=t(k,i)*sine+t(k,j)*cosine
t(k,j)=q t(k,i)=p
40 continue t(k,j)=q
45 d=dabs(sine) enddo
if (d.le.amax) goto 50 45 d=dabs(sine)
imax=i if (d.le.amax) goto 50
jmax=j imax=i
amax=d jmax=j
dmax=delta amax=d
50 continue dmax=delta
60 continue 50 continue
end do
end do
70 format (' iter=',i4,' largest rotation=',f12.8, 70 format (' iter=',i4,' largest rotation=',f12.8,
* ', vectors',i3,' and',i3,', incr. of diag. squares=',g12.5) * ', vectors',i3,' and',i3,', incr. of diag. squares=',g12.5)
71 format (' i,j,a,b,sin,cos,delta =',2i3,5f10.5) 71 format (' i,j,a,b,sin,cos,delta =',2i3,5f10.5)
@ -285,22 +301,27 @@ C conv=1.d-6
* 'in subroutine maxovl ***'//) * 'in subroutine maxovl ***'//)
stop stop
100 continue 100 continue
do 120 j=1,n do j=1,n
if (s(j,j).gt.0.d0) goto 120 if (s(j,j).gt.0.d0) cycle
do 105 i=1,m do i=1,m
105 s(i,j)=-s(i,j) s(i,j)=-s(i,j)
do 110 i=1,n enddo
110 t(i,j)=-t(i,j) do i=1,n
120 continue t(i,j)=-t(i,j)
enddo
enddo
sum=0.d0 sum=0.d0
do 125 i=1,m do i=1,m
125 sum=sum+s(i,i)*s(i,i) sum=sum+s(i,i)*s(i,i)
enddo
sum=sum/m sum=sum/m
do 122 i=1,m do i=1,m
do 122 j=1,n do j=1,n
sw=s(i,j) sw=s(i,j)
s(i,j)=w(i,j) s(i,j)=w(i,j)
122 w(i,j)=sw w(i,j)=sw
enddo
enddo
return return
end end