yapanuwan
/
BergEC-jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Remove OSBRACKETS from repo

This commit is contained in:
ysyapa 2024-04-16 19:04:49 +00:00
parent 0ff9786bef
commit 3cc2b39d41
5 changed files with 0 additions and 1515 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

71
OSBRACKETS/README
View File

@ -1,71 +0,0 @@
The programs are written in Fortran-90. The double precision is adopted.
The following six files make up the package:
allosbrac.f90
testallosbrac.f90
allosoutput
osbrac.f90
testosbrac.f90
osoutput
These are two independent groups of files that correspond to the two versions of the present
program.
In one version, the subroutine producing the brackets is named ALLOSBRAC. It is contained in
the file allosbrac.f90 along with the routines this subroutine uses.
In the other version, the subroutine producing the brackets is named OSBRAC. It is contained
in the file osbrac.f90 along with the routines this subroutine uses.
To test the set of programs contained in the file allosbrac.f90, one may say, e.g.,
gfortran -O2 testallosbrac.f90 allosbrac.f90. To test the set of programs contained in the
file osbrac.f90, one may say, e.g., gfortran -O2 testosbrac.f90 osbrac.f90. The files
allosoutput and osoutput are respective output files that contain results of the tests.
The files testallosbrac.f90 and testosbrac.f90 contain, respectively, only the programs
TESTALLOSBRAC and TESTOSBRAC. The tests performed are described in comment lines
in these programs and in the accompanying CPC paper.
The parameters that are set at the moment in the programs TESTALLOSBRAC and TESTOSBRAC
are just those for which the results in the files allosoutput and osoutput are listed. With these
parameters, the programs run less than a second on a notebook.
-------------------------------------------------------------------------------------------------------------------------------------------
As said above, the brackets are produced either by the subroutine ALLOSBRAC or by the
subroutine OSBRAC. The subroutines contain generous comments.
-------------------------------------------------------------------------------------------------------------------------------------------
In the case when ALLOSBRAC is employed, all the routines used to calculate the brackets,
which are contained in the file allosbrac.f90, are as follows: the subroutine ALLOSBRAC (with
comments), the subroutines ARR, FLPHI (with comments), COEFREL (with comments), and the
function WIGMOD.
The subroutine ALLOSBRAC calls for the ARR, COE, FLPHI, and COEFREL subroutines. FLPHI calls
for the function WIGMOD. The meaning of the parameters of the ALLOSBRAC subroutine is
explained in comments at its beginning.
------------------------------------------------------------------------------------------------------------------------------------------
In the case when OSBRAC is employed, all the routines used to calculate the brackets,
which are contained in the file osbrac.f90, are as follows: the subroutine OSBRAC (with
comments), the subroutines ARR, FLPHI (with comments), COEFREL (with comments), and the
function WIGMOD.
The subroutine OSBRAC calls for the ARR, COE, FLPHI, and COEFREL subroutines. FLPHI calls
for the function WIGMOD. The meaning of the parameters of the OSBRAC subroutine is explained
in comments at its beginning.
-------------------------------------------------------------------------------------------------------------------------------------------
While it is clear from the comments in the beginnings of ALLOSBRAC and OSBRAC how to
implement these subroutines, this is also obvious from examples listed in the above mentioned
TESTALLOSBRAC and TESTOSBRAC programs.
In general, to compile codes that include the present programs one may say e.g.,
f95 main.f90 ... allosbrac.f90 ... or, alternatively, f95 main.f90 ... osbrac.f90 ....

451
OSBRACKETS/allosbrac.f90
View File

@ -1,451 +0,0 @@
SUBROUTINE ALLOSBRAC(NQMAX,LMIN,LMAX,CO,SI,BRAC)
! COMMENTS ON USE:
!---------------------------------------------------------------------------------------------------------------------------------------
! THE OSCILLATOR BRACKETS <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\VARPHI ARE
! CALCULATED HERE. THIS IS DONE WITH THE HELP OF EQS. (18), (20)-TYPE, AND (22) IN
! THE ACCOMPANYING CPC PAPER. THE QUANTITIES L1, L2, L1P, AND L2P IN THE ABOVE
! DEFINITION OF THE BRACKETS ARE THE PARTIAL ANGULAR MOMENTA, L IS THE TOTAL
! ANGULAR MOMENTUM, AND N1, N2, N1P, N2P ARE THE RADIAL QUANTUM NUMBERS.
! IN THE NOTATION LIKE L1P, ETC., "P" SYMBOLIZES "PRIMED" HERE AND BELOW.
! THE SUBROUTINE RETURNS THE ARRAY OF ALL THE BRACKETS SUCH THAT
! L1+L2+2*(N1+N2).LE.NQMAX, LMIN \LE L \LE LMAX, AND THE BRACKETS PERTAIN TO STATES
! OF THE SAME PARITY (-1)^(L1+L2) WHICH IS THE PARITY OF NQMAX.
! THE L1, L2, L1P, AND L2P ORBITAL MOMENTA ARE EXPRESSED IN THE SUBROUTINE IN
! TERMS OF THE M, N, MP, AND NP VARIABLES DEFINED AS FOLLOWS
! M=(L1+L2-L-NN)/2, N=(L1-L2+L-NN)/2, MP=(L1P+L2P-L-NN)/2, NP=(L1P-L2P+L-NN)/2
! WHERE NN EQUALS 0 OR 1 WHEN, RESPECTIVELY, NQMAX-L IS EVEN OR ODD. ONE THEN HAS
! L1 = M+N+NN, L2 = M-N+L, L1P = MP+NP+NN, L2P = MP-NP+L.
! WHEN L1, L2, L1P AND L2P TAKE ALL THE VALUES ALLOWED AT GIVEN L AND GIVEN PARITY
! THE N AND NP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO L-NN AND
! THE M AND MP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO (NQMAX-L-NN)/2.
! ONE ALSO HAS NQ=2*(M+N1+N2)+L+NN=2*(MP+N1P+N2P)+L+NN.
! ALL THE PARAMETERS OF THE SUBROUTINE BUT BRAC ARE INPUT ONES. THE MEANING
! OF THE NQMAX, LMIN, AND LMAX PARAMETERS IS SEEN FROM THE ABOVE DESCRIPTION.
! THE BRAC PARAMETER DESIGNATES THE ARRAY OF OUTPUT BRACKETS. IT IS OF THE FORM
! BRAC(NP,N1P,MP,N1,N2,N,M,L) . (AS SAID ABOVE, L1=L1(N,M), L2=L2(N,M), L1P=L1P(NP.MP),
! AND L2P=L2P(NP,MP). THE QUANTITY N2P IS DETERMINED BY THE EQUALITY
! MP+N1P+N2P=M+N1+N2.) THE ORDER OF THE ARGUMENTS OF BRAC CORRESPONDS TO
! NESTING OF THE LOOPS AT ITS COMPUTATION.
! THE ROUTINE PARAMETERS CO AND SI ARE AS FOLLOWS, CO = COS(PHI) AND SI = SIN(PHI).
! THESE QUANTITIES DEFINE THE PSEUDO(!!!)ORTHOGONAL TRANSFORMATION
! XI1P = CO*XI1+SI*XI2, XI2P = SI*XI1-CO*XI2.
! (BRACKETS PERTAINING TO THE CASE OF THE ORTHOGONAL TRANSFORMATION
! XI1P = CO*XI1+SI*XI2, XI2P = -SI*XI1+CO*XI2 ARE SIMPLY EXPRESSED IN TERMS OF
! THOSE CALCULATED IN THE PRESENT PROGRAM, SEE THE ACCOMPANYING CPC PAPER.
! THE MEANING OF THE ARRAYS, OTHER THAN BRAC, ENTERING THE DIMENSION LIST IS AS
! FOLLOWS. B IS A SUBSIDIARY ARRAY TO PERFORM THE RECURSION.
! FAC(I)=I!, DFAC(I)=(2I+1)!!, AND DEFAC(I)=(2I)!!. RFAC(I)=SQRT((2I)!)/I!
! THESE ARRAYS ARE PRODUCED BY THE "ARR" ROUTINE WHICH IS CALLED FROM THE PRESENT
! ROUTINE AND WHICH IS CONTAINED IN THE PRESENT FILE. THE SET UPPER BOUNDS OF THESE
! ARRAYS ARE SUFFICIENT FOR THE COMPUTATION. THESE ARRAYS ARE USED IN THE PRESENT
! ROUTINE AND IN THE "FLPHI" ROUTINE.
! THE A ARRAY IS A(N,L)=(-1)**N/SQRT(DEFAC(N)*DFAC(N+L)). IT APPEARS BOTH IN THE N1=N2=0
! BRACKETS AND IN THE RELATION BETWEEN THE < | > AND [ | ] TYPE BRACKETS.
! THE BI ARRAY IS BI(M,N)=FAC(N)/[FAC(M)*FAC(N-M)]. IT ENTERS THE 3J SYMBOLS. IT IS USED
! IN THE "FLPHI" ROUTINE AND IN THE FUNCTION WIGMOD.
! THE FL(NP,MP,N) ARRAY REPRESENTS THE QUANTITY
! [(2*L1P+1)*(2*L2P+1)]^{1/2}*F_L^VARPHI WHERE F_L^VARPHI IS GIVEN BY EQ. (23) IN THE
! ACCOMPANYING CPC PAPER. THIS ARRAY IS PRODUCED IN ADVANCE BY THE
! "FLPHI" ROUTINE WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS
! CONTAINED IN THE PRESENT FILE. THE VARIABLES MP, NP, AND N ARE DEFINED ABOVE ALONG
! WITH THEIR UPPER BOUNDS.
! THE PSIP(P,Q) AND PSIM(P,Q) ARRAYS REPRESENT THE QUANTITIES (4) IN THE ACCOMPANYING
! CPC PAPER. THEY ARE PRODUCED IN ADVANCE BY THE "COEFREC" ROUTINE
! WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS CONTAINED IN THE PRESENT
! FILE. THEIR ARGUMENTS P AND Q TAKE THE VALUES L1P, L1P+/-1 AND L2P, L2P+/-1,
! RESPECTIVELY. THE BOUNDS OF THE PSIP AND PSIM ARRAYS ARE SUCH THAT ALL THE P AND Q
! VALUES REQUIRED TO PERFORM THE RECURSION ARE PROVIDED.
! THE MENTIONED "FLPHI" ROUTINE USES THE FUNCTION WIGMOD WHICH IS ALSO CONTAINED
! IN THE PRESENT FILE.
!----------------------------------------------------------------------------------------------------------------------------------------
DOUBLE PRECISION BRAC,FL,PSIP,PSIM,A,FAC,DFAC,DEFAC,RFAC,BI,CO,SI,T,CO2,SI2,SC,FA,&
PLFACT,D2L12P,TN,PA,S,FA1,FA2
DIMENSION&
BRAC(0:LMAX,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,&
0:(NQMAX-LMIN)/2,0:LMAX,0:(NQMAX-LMIN)/2,LMIN:LMAX),&
! THE ARRAY IS OF THE FORM BRAC(NP,N1P,MP,N1,N2,N,M,L)
FL(0:LMAX,0:(NQMAX-LMIN)/2,0:LMAX),&
PSIP((NQMAX+LMAX)/2+1,(NQMAX+LMAX)/2+1),&
PSIM((NQMAX+LMAX)/2+1,(NQMAX+LMAX)/2+1),A(0:NQMAX/2+1,0:NQMAX),&
FAC(0:2*NQMAX+1),DFAC(0:NQMAX+1),DEFAC(0:NQMAX/2+1),RFAC(0:NQMAX),&
BI(0:2*NQMAX+1,0:2*NQMAX+1)
IF (NQMAX.GT.84) THEN
WRITE(6,*)'IN THE R(*8) COMPUTATION NQMAX SHOULD NOT EXCEED 84'
STOP
ENDIF
IF (CO.EQ.0.D0.OR.SI.EQ.0.D0) THEN
WRITE(6,*)'THE PROGRAM IS OF NO USE AT ZERO COS(PHI) OR SIN(PHI) VALUES'
WRITE(6,*)'COS(PHI)=',CO,' SIN(PHI)=',SI
STOP
ENDIF
IF (NQMAX.LT.LMAX) THEN
WRITE(6,*)'NQMAX=',NQMAX,' L=',LMAX
WRITE(6,*)'L SHOULD NOT EXCEED NQMAX'
STOP
ENDIF
T=SI/CO
CO2=CO**2
SI2=SI**2
SC=SI*CO
CALL ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
CALL COE(2*NQMAX+1,FAC,BI)
! THE A_{NL} ARRAY:
DO NI=0,NQMAX/2+1
NA=NI-2*(NI/2)
K=1
IF(NA.EQ.1)K=-1
DO LI=0,NQMAX-NI
A(NI,LI)=K/SQRT(DEFAC(NI)*DFAC(NI+LI))
ENDDO
ENDDO
! BRACKETS [N_1'L_1'N_2'L_2'|0L_10L_2]_L^\PHI TO START THE RECURSION, EQ. (22) IN
! THE ACCOMPANYING CPC PAPER.
! THESE BRACKETS ARE REPRESENTED AS BRAC(NP,N1P,MP,0,0,N,M,L).
! SUBSIDIARY QUANTITIES:
DO L=LMIN,LMAX
CALL FLPHI(NQMAX,LMAX,LMIN,L,T,BI,RFAC,FL)
! THIS SUBROUTINE PRODUCES THE ARRAY FL USED BELOW.
NQML=NQMAX-L
NQMLD2=NQML/2
NN=NQML-2*(NQMLD2)
CALL COEFREC((NQMAX+LMAX)/2+1,L,NN,NQMAX,PSIP,PSIM)
NMAX=L-NN
MMAX=(NQML-NN)/2
DO M=0,MMAX
DO N=0,NMAX
L1=M+N+NN
L2=M-N+L
L1L2=L1+L2
L1L2ML=L1L2-L
FA=SQRT((2.D0*L1+1)*(2*L2+1))*CO**L1L2
PLFACT=SQRT(FAC(L1L2+L+1)*FAC(L1L2ML))
! COMPUTATION OF THE N1=N2=0 BRACKETS. THEY ARE REPRESENTED AS
! BRAC(NP,N1P,MP,0,0,N,M,L).
MPMAX=(L1L2ML-NN)/2
DO MP=0,MPMAX
L12P=2*MP+L+NN
! THIS IS BY DEFINITION, L12P=L1P+L2P
D2L12P=1.D0/2.D0**L12P
N12P=MPMAX-MP
! INDEED, 2*N12P=2*(N1P+N2P)=NQ-L12P. FOR THE BRACKETS WE COMPUTE NOW
! WE HAVE NQ=L1L2. THEN 2*N12P=L1L2-L12P. SINCE L1L2=2*MPMAX-L-NN AND
! L12P=2*MP+L+NN, ONE GETS N12P=MPMAX-MP.
DO N1P=0,N12P
N1PA=N1P-2*(N1P/2)
IF (N1PA.EQ.0) THEN
K=1
ELSE
K=-1
ENDIF
! K=(-1)**N1P
N2P=N12P-N1P
TN=T**N12P*K
DO NP=0,NMAX
L1P=MP+NP+NN
L2P=L12P-L1P
! BY DEFINITION ONE ALSO HAS: L2P=L2P(MP,NP)=MP-NP+L
PA=A(N1P,L1P)*A(N2P,L2P)
BRAC(NP,N1P,MP,0,0,N,M,L)=FL(NP,MP,N)*FA*TN*D2L12P&
*PLFACT*PA*PA
ENDDO ! NP
ENDDO ! N1P
ENDDO ! MP
! RECURSION TO OBTAIN THE GENERAL FORM BRACKETS
! [N_1'L_1'N_2'L_2'|N_1L_1N_2L_2]_L^\PHI.
N12MAX=(NQMAX-L1L2)/2
MPMAX0=MPMAX
! THE N1=0, N2-1-->N2 RECURSION, EQ. (20) IN THE ACCOMPANYING CPC PAPER WITH THE
! MODIFICATION POINTED OUT THERE:
DO N2=1,N12MAX
MPMAX=MPMAX+1
! THIS MPMAX VALUE IS (NQ-L-NN)/2 WHERE NQ=L1+L2+2*N2.
DO MP=0,MPMAX
N12P=MPMAX-MP
DO N1P=0,N12P
! MP=(L1P+L2P-L-NN)/2. THEREFORE, N12P=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.
DO NP=0,NMAX
L1P=MP+NP+NN
L2P=MP-NP+L
S=0.D0
! BELOW THE RESTRICTIONS ARE IMPOSED ON THE BRAC ARRAY ENTERING THE RIGHT-HAND
! SIDE OF THE RECURRENCE FORMULAE. THIS ARRAY IS OF THE FORM BRAC(K1,K2,K3,...) WHERE
! K1=K1(NP), K2=N2(N1P), AND K3=K3(MP). NAMELY, K1=NP, OR NP+1, OR NP-1;
! K2=N1P OR N1P-1; K3=MP, OR MP-1, OR MP+1. (THE COMBINATION K2=N1P AND K3=MP+1
! DOES NOT ARISE IN THE RECURRENCE FORMULAE.) THE QUANTITIES K1, K2, AND K3
! REPRESENT, RESPECTIELY, NP, N1P, AND MP VALUES AT THE PRECEDING STAGE OF THE
! RECURSION. THE RESTRICTIONS ENSURE THAT K1, K2, AND K3 RANGE WITHIN THE LIMITS
! PERTAINING TO THE BRAC ARRAY OBTAINED AT THAT PRECEDING STAGE OF THE RECURSION.
! THUS K1, K2, AND K3 SHOULD BE NON-NEGATIVE. THEREFORE, WHEN K1=NP-1 THE VALUE
! NP=0 IS TO BE EXCLUDED, WHEN K2=N1P-1 THE VALUE N1P=0 IS TO BE EXCLUDED, AND
! WHEN K3=MP-1 THE VALUE MP=0 IS TO BE EXCLUDED.
! FURTHERMORE, IT SHOULD BE K1 \LE NMAX AND AT THE SAME TIME ONE HAS NP \LE NMAX.
! THEREFORE, WHEN K1=NP+1 THE VALUE NP=NPMAX IS TO BE EXCLUDED.
! IN ADDITION, IT SHOULD BE K2+K3 \LE MPMAX-1 WHILE ONE HAS N1P+MP \LE MPMAX. WHEN
! K2=N1P AND K3=MP-1, OR K2=N1P-1 AND K3=MP-1, OR K2=N1P-1 AND K3=MP ONE
! AUTOMATICALLY HAS N1P+MP \LE MPMAX. THEREFORE, IN THESE CASES THE CONDITION
! K2+K3 \LE MPMAX-1 DOES NOT CREATE ANY RESTRICTIONS.
! BUT WHEN K2=N1P AND K3=MP, OR K2=N1P-1 AND K3=MP+1 THE CASE N1P+MP=MPMAX
! IS TO BE FORBIDDEN.
! IN THE RECURRENCE FORMULAE BELOW THE DESCRIBED RESTRICTIONS ARE
! IMPOSED.
IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,0,N2-1,N,M,L)*PSIP(L1P,L2P)
IF (N1P.NE.0.AND.N1P.NE.N12P) S=S+&
BRAC(NP,N1P-1,MP+1,0,N2-1,N,M,L)*PSIP(L1P+1,L2P+1)
! RECALL THAT N12P=MPMAX-MP.
IF (NP.NE.NMAX.AND.N1P.NE.0) S=S-&
BRAC(NP+1,N1P-1,MP,0,N2-1,N,M,L)*PSIM(L1P+1,L2P)
IF (NP.NE.0.AND.N1P.NE.N12P) S=S-&
BRAC(NP-1,N1P,MP,0,N2-1,N,M,L)*PSIM(L1P,L2P+1)
S=S*SC
IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,0,N2-1,N,M,L)*SI2
IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,0,N2-1,N,M,L)*CO2
BRAC(NP,N1P,MP,0,N2,N,M,L)=S
ENDDO ! NP
ENDDO ! N1P
ENDDO ! MP
ENDDO ! N2
! N1-1-->N1 RECURSION, EQ. (20) IN THE ACCOMPANYING CPC PAPER:
DO N2=0,N12MAX
MPMAX=MPMAX0+N2
! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2=MPMAX0+N2 SINCE NQ=L1+L2+2*N2.
! THE RECURSION:
DO N1=1,N12MAX-N2
MPMAX=MPMAX+1
! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2 AND NQ=L1+L2+2*N2+2*N1.
DO MP=0,MPMAX
N12P=MPMAX-MP
DO N1P=0,N12P
! MP=(L1P+L2P-L-NN)/2. THEREFORE, N1PMAX=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.
DO NP=0,NMAX
L1P=MP+NP+NN
L2P=MP-NP+L
S=0.D0
IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,N1-1,N2,N,M,L)*PSIP(L1P,L2P)
IF (N1P.NE.0.AND.N1P.NE.N12P) S=S+&
BRAC(NP,N1P-1,MP+1,N1-1,N2,N,M,L)*PSIP(L1P+1,L2P+1)
IF (NP.NE.NMAX.AND.N1P.NE.0) S=S-&
BRAC(NP+1,N1P-1,MP,N1-1,N2,N,M,L)*PSIM(L1P+1,L2P)
IF (NP.NE.0.AND.N1P.NE.N12P) S=S-&
BRAC(NP-1,N1P,MP,N1-1,N2,N,M,L)*PSIM(L1P,L2P+1)
! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
S=-S*SC
IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,N1-1,N2,N,M,L)*CO2
IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,N1-1,N2,N,M,L)*SI2
! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
BRAC(NP,N1P,MP,N1,N2,N,M,L)=S
ENDDO ! NP
ENDDO ! N1P
ENDDO ! MP
ENDDO ! N1
ENDDO ! N2
! RENORMALIZATION OF THE BRACKETS, EQ. (18) IN THE ACCOMPANYING CPC PAPER:
DO N2=0,N12MAX
MPMAX1=MPMAX0+N2
DO N1=0,N12MAX-N2
MPMAX=MPMAX1+N1
FA1=A(N1,L1)*A(N2,L2)
DO MP=0,MPMAX
N12P=MPMAX-MP
DO N1P=0,MPMAX-MP
N2P=N12P-N1P
DO NP=0,NMAX
L1P=MP+NP+NN
L2P=MP-NP+L
FA2=A(N1P,L1P)*A(N2P,L2P)
BRAC(NP,N1P,MP,N1,N2,N,M,L)=BRAC(NP,N1P,MP,N1,N2,N,M,L) &
*FA1/FA2
ENDDO ! NP
ENDDO ! N1P
ENDDO ! MP
ENDDO ! N1
ENDDO ! N2
ENDDO ! N
ENDDO ! M
ENDDO ! L
RETURN
END
SUBROUTINE ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
! FAC(I), DFAC(I),DEFAC(I), AND RFAC(I) ARE, RESPECTIVELY, THE QUANTITIES
! I!, (2I+1)!!, (2I)!!, AND SQRT((2*I)!)/I!
DOUBLE PRECISION FAC,DFAC,DEFAC,RFAC
DIMENSION FAC(0:2*NQMAX+1),DFAC(0:NQMAX+1),DEFAC(0:NQMAX/2+1),RFAC(0:NQMAX)
FAC(0)=1.D0
DFAC(0)=1.D0
DEFAC(0)=1.D0
RFAC(0)=1.D0
DO I=1,2*NQMAX+1
FAC(I)=FAC(I-1)*I
ENDDO
DO I=1,NQMAX+1
DFAC(I)=DFAC(I-1)*(2*I+1)
ENDDO
DO I=1,NQMAX/2+1
DEFAC(I)=DEFAC(I-1)*2*I
ENDDO
DO I=1,NQMAX
RFAC(I)=RFAC(I-1)*2*SQRT(1-0.5D0/I)
ENDDO
RETURN
END
SUBROUTINE FLPHI(NQMAX,LMAX,LMIN,L,T,BI,RFAC,FL)
! PROVIDES THE QUANTITY SQRT((2L1P+1)*(2L2P+1))*F_L^\VARPHI, SEE EQ. (23) IN THE
! ACCOMPANYING PAPER, IN THE FORM OF THE FL ARRAY.
! USES THE FUNCTION WIGMOD.
! THIS FLPHI DIFFERS FROM THAT IN THE OTHER FILE osbrac.f90.
DOUBLE PRECISION BI,RFAC,FL,T,T2,SQP,F,WIGMOD,z
DIMENSION&
FL(0:LMAX,0:(NQMAX-LMIN)/2,0:LMAX),BI(0:2*NQMAX+1,0:2*NQMAX+1),RFAC(0:NQMAX)
! THE OUTPUT IS FL(NP,MP,N) WHERE MP=(L1P+L2P-L-NN)/2, L1P+L2P=L1T+L2T,
! NP=(L1P-L2P+L-NN)/2, AND N=(L1-L2+L-NN)/2, L1-L2=L1T-L2T.
T2=T*T
NQMAL=NQMAX-L
MPMAX=NQMAL/2
NN=NQMAL-2*MPMAX
LMNN=L-NN
LL3=2*L
DO N=0,LMNN
L1ML2=2*N-LMNN
DO MP=0,MPMAX
L1PL2P=2*MP+L+NN
DO NP=0,LMNN
L1PML2P=2*NP-LMNN
L1P=(L1PL2P+L1PML2P)/2
L2P=(L1PL2P-L1PML2P)/2
M3=L2P-L1P
SQP=SQRT((2*L1P+1.D0)*(2*L2P+1))
L1T=(L1PL2P+L1ML2)/2
L2T=(L1PL2P-L1ML2)/2
N3=L1T+L2T-L
IMIN=ABS(L1T-L1P)
IMAX=MIN(L1P+L1T,L2P+L2T)
NAL1=(L1T+L1P-IMIN)/2
NAL2=(L1T-L1P+IMIN)/2
NAL3=(L2T-L2P+IMIN)/2
NAL4=(L2T+L2P-IMIN)/2
NAL4P=NAL4-2*(NAL4/2)
F=0.D0
z=t**imin
if(nal4p.ne.0)z=-z
DO I=IMIN,IMAX,2
F=F+RFAC(NAL1)*RFAC(NAL2)*RFAC(NAL3)*RFAC(NAL4)*&
WIGMOD(L1T,L2T,L,L1P-I,I-L2P,BI,2*NQMAX+1)*z
NAL1=NAL1-1
NAL2=NAL2+1
NAL3=NAL3+1
NAL4=NAL4-1
z=-z*t2
ENDDO
FL(NP,MP,N)=F*SQRT((2*L1P+1)*(2*L2P+1)&
*BI(L1T+L-L2T,2*L1T)*BI(N3,2*L2T)/((LL3+1)*BI(N3,L1T+L2T+L+1)*BI(L+M3,LL3)))
ENDDO
ENDDO
ENDDO
RETURN
END
FUNCTION WIGMOD(L1,L2,L3,M1,M2,BI,NMAX)
! THE 3J SYMBOL IN TERMS OF BINOMIAL COEFFICIENTS WITHOUT THE FACTOR
! SQRT(F) WHERE F IS AS FOLLOWS,
! F=BI(L1+L3-L2,2*L1)*BI(N3,2*L2)/((LL3+1)*BI(N3,L1+L2+L3+1)*BI(L3+M3,LL3))
! WITH N3=L1+L2-L3, LL3=2*L3, AND M3=-M1-M2.
DOUBLE PRECISION WIGMOD,BI,S
DIMENSION BI(0:NMAX,0:NMAX)
M3=-M1-M2
N3=L1+L2-L3
LM1=L1-M1
LP2=L2+M2
KMIN=MAX(0,L1-L3+M2,L2-L3-M1)
KMAX=MIN(N3,LM1,LP2)
S=0.d0
NPH=1
DO K=KMIN,KMAX
S=S+NPH*BI(K,N3)*BI(LM1-K,L1+L3-L2)*BI(LP2-K,L2+L3-L1)
NPH=-NPH
ENDDO
WIGMOD=S/SQRT(BI(LM1,2*L1)*BI(LP2,2*L2))
NY=KMIN+L1-L2-M3
NYP=NY-2*(NY/2)
IF(NYP.NE.0)WIGMOD=-WIGMOD
RETURN
END
SUBROUTINE COE(NMAX,FAC,BI)
DOUBLE PRECISION FAC,BI
DIMENSION FAC(0:NMAX),BI(0:NMAX,0:NMAX)
DO N=0,NMAX
DO M=0,N/2
BI(M,N)=FAC(N)/(FAC(M)*FAC(N-M))
BI(N-M,N)=BI(M,N)
ENDDO
ENDDO
RETURN
END
SUBROUTINE COEFREC(MAXCOE,L,NN,NQMAX,PSIP,PSIM)
! CALCULATES THE COEFFICIENTS OF THE RECURSION FORMULA, EQ. (21) IN THE ACCOMPANYING
! CPC PAPER
! MAXCOE EQUALS INT((NQMAX+LMAX)/2)+1 AT CALLS OF THIS ROUTINE.
DOUBLE PRECISION PSIP,PSIM
DIMENSION PSIP(MAXCOE,MAXCOE),PSIM(MAXCOE,MAXCOE)
LP=L+1
LM=L-1
DO M1=1,(NQMAX+L)/2+1
DO M2=1,M1
M1P2=M1+M2
NA=M1P2+L
NNA=NA-2*(NA/2)
M1M2=M1-M2
! CALCULATION OF PSIP. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P AND L2P OR
! L1P+1 AND L2P+1.)
IF (NNA.EQ.NN.AND.M1P2.GE.L.AND.M1P2.LE.NQMAX.AND.&
ABS(M1M2).LE.L) THEN
IF (M1P2.GT.LP) THEN
M1P2L=M1P2+L
M1P2ML=M1P2-L
PSIP(M1,M2)=SQRT(M1P2L*(M1P2L+1.D0)*M1P2ML*(M1P2ML-1)/((4*M1*M1-1)&
*(4*M2*M2-1)))
ELSE
PSIP(M1,M2)=0.D0
ENDIF
PSIP(M2,M1)=PSIP(M1,M2)
ENDIF
! CALCULATION OF PSIM. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P+1 AND L2P
! OR L2P+1 AND L1P.)
IF (NNA.NE.NN.AND.M1P2.GE.LP.AND.M1P2.LE.NQMAX&
.AND.M1M2.LE.LP.AND.M1M2.GE.-LM) THEN
IF (M1M2.LT.L) THEN
M1M2L=M1M2+L
M1M2ML=M1M2-L
PSIM(M1,M2)=SQRT(M1M2L*(M1M2L+1.D0)*M1M2ML*(M1M2ML-1)/&
((4*M1*M1-1)*(4*M2*M2-1)))
ELSE
PSIM(M1,M2)=0.D0
ENDIF
PSIM(M2,M1)=PSIM(M1,M2)
ENDIF
ENDDO
ENDDO
RETURN
END

467
OSBRACKETS/osbrac.f90
View File

@ -1,467 +0,0 @@
SUBROUTINE OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
! COMMENTS ON USE:
!---------------------------------------------------------------------------------------------------------------------------------------
! THE OSCILLATOR BRACKETS <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\VARPHI ARE
! CALCULATED HERE. THIS IS DONE WITH THE HELP OF EQS. (18), (20)-TYPE, AND (22) IN
! THE ACCOMPANYING CPC PAPER. THE QUANTITIES L1, L2, L1P, AND L2P IN THE ABOVE
! DEFINITION OF THE BRACKETS ARE THE PARTIAL ANGULAR MOMENTA, L IS THE TOTAL
! ANGULAR MOMENTUM, AND N1, N2, N1P, N2P ARE THE RADIAL QUANTUM NUMBERS.
! IN THE NOTATION LIKE L1P, ETC., "P" SYMBOLIZES "PRIMED" HERE AND BELOW.
! N,M,L, AND NQMAX ARE INPUT PARAMETERS OF THE SUBROUTINE, AND BRAC IS THE
! OUTPUT PARAMETER REPRESENTING THE ARRAY OF CALCULATED BRACKETS. THE PARAMETERS
! N AND M REPRESENT THE L1 AND L2 QUANTUM NUMBERS, SEE BELOW, SO THAT L1=L1(N,M)
! AND L2=L2(N,M).
! LET US USE THE NOTATION NQ=L1+L2+2*(N1+N2). THE ARRAY BRAC CONTAINS ALL BRACKETS
! EXISTING AT GIVEN INPUT L, L1, AND L2 VALUES AND HAVING N1 AND N2 SUCH THAT
! NQ.LE.NQMAX. IT IS CLEAR THAT NQ AND NQMAX PERTAINING TO ALL THESE BRACKETS ARE OF
! THE SAME PARITY EQUAL TO (-1)^(L1+L2).
! IN OTHER WORDS, THE COMPUTED BRACKETS PERTAIN TO STATES WITH GIVEN N, M, AND L
! VALUES SUCH THAT NQ.LE.NQMAX. AUTOMATICALLY, THESE STATES ARE OF THE SAME PARITY
! AS THE PARITY OF NQMAX. (AND THIS PARITY EQUALS (-1)^(L1+L2).)
! THE L1, L2, L1P, AND L2P ORBITAL MOMENTA ARE EXPRESSED IN TERMS OF THE N, M, NP, AND
! MP VARIABLES AS FOLLOWS
! N=(L1-L2+L-NN)/2, M=(L1+L2-L-NN)/2, NP=(L1P-L2P+L-NN)/2, MP=(L1P+L2P-L-NN)/2
! WHERE NN EQUALS 0 OR 1 WHEN, RESPECTIVELY, NQMAX-L IS EVEN OR ODD. ONE THEN HAS
! L1 = M+N+NN, L2 = M-N+L, L1P = MP+NP+NN, L2P = MP-NP+L.
! WHEN L1, L2, L1P AND L2P TAKE ALL THE VALUES ALLOWED AT GIVEN L AND GIVEN PARITY
! THE N AND NP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO L-NN AND
! THE M AND MP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO (NQMAX-L-NN)/2.
! ONE ALSO HAS NQ=2*(M+N1+N2)+L+NN=2*(MP+N1P+N2P)+L+NN.
! THE SAID ABOVE DETERMINES THE MEANING OF THE N,M, AND NQMAX PARAMETERS OF
! THE SUBROUTINE.
! THE BRAC ARRAY IS OF THE FORM BRAC(NP,N1P,MP,N1,N2). AS SAID ABOVE, NP AND MP
! RERESENT L1P AND L2P. THE MISSING N2P VALUE IS DETERMNED BY THE EQUALITY
! MP+N1P+N2P=M+N1+N2 THAT IS ANOTHER FORM OF THE RELATION
! 2(N1+N2)+L1+L2=2(N1P+N2P)+L1P+L2P.
! THE ROUTINE PARAMETERS CO AND SI ARE AS FOLLOWS, CO = COS(PHI) AND SI = SIN(PHI).
! THESE QUANTITIES DEFINE THE PSEUDO(!)ORTHOGONAL TRANSFORMATION
! XI1P = CO*XI1+SI*XI2, XI2P = SI*XI1-CO*XI2.
! (BRACKETS PERTAINING TO THE CASE OF THE ORTHOGONAL TRANSFORMATION
! XI1P = CO*XI1+SI*XI2, XI2P = -SI*XI1+CO*XI2 ARE SIMPLY EXPRESSED IN TERMS OF
! THOSE CALCULATED IN THE PRESENT PROGRAM, SEE THE ACCOMPANYING CPC PAPER.
! THE MEANING OF THE FIRSTCALL PARAMETER IS EXPLAINED BELOW.
! THE MEANING OF THE ARRAYS, OTHER THAN BRAC, ENTERING THE DIMENSION LIST IS AS
! FOLLOWS. B IS A SUBSIDIARY ARRAY TO PERFORM THE RECURSION.
! FAC(I)=I!, DFAC(I)=(2I+1)!!, AND DEFAC(I)=(2I)!!, AND RFAC(I)=SQRT((2I)!)/I!
! IS A SUBSIDIARY ARRAY TO CALCULATE THE FL ARRAY, SEE BELOW. THESE ARRAYS ARE
! PRODUCED BY THE "ARR" ROUTINE WHICH IS CALLED FROM THE PRESENT ROUTINE AND
! WHICH IS CONTAINED IN THE PRESENT FILE. THE SET UPPER BOUNDS OF THESE
! ARRAYS ARE SUFFICIENT FOR THE COMPUTATION. THESE ARRAYS ARE USED IN THE PRESENT
! ROUTINE AND IN THE "FLPHI" ROUTINE.
! THE A ARRAY IS A(N,L)=(-1)**N/SQRT(DEFAC(N)*DFAC(N+L)). IT APPEARS BOTH
! IN THE N1=N2=0 BRACKETS AND IN THE RELATION BETWEEN THE < | > AND [ | ] TYPE
! BRACKETS.
! THE BI ARRAY IS BI(M,N)=FAC(N)/[FAC(M)*FAC(N-M)]. IT ENTERS THE 3J SYMBOLS. IT IS USED
! IN THE "FLPHI" ROUTINE AND IN THE FUNCTION WIGMOD.
! THE FL(NP,MP,N) ARRAY REPRESENTS THE QUANTITY
! [(2*L1P+1)*(2*L2P+1)]^{1/2}*F_L^VARPHI WHERE F_L^VARPHI IS GIVEN BY EQ. (23) IN THE
! ACCOMPANYING CPC PAPER. THIS ARRAY IS PRODUCED IN ADVANCE BY THE
! "FLPHI" ROUTINE WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS
! CONTAINED IN THE PRESENT FILE. THE VARIABLES MP, NP, AND N ARE DEFINED ABOVE ALONG
! WITH THEIR UPPER BOUNDS.
! THE PSIP(P,Q) AND PSIM(P,Q) ARRAYS REPRESENT THE QUANTITIES (4) IN THE ACCOMPANYING
! CPC PAPER. THEY ARE PRODUCED IN ADVANCE BY THE "COEFREC" ROUTINE
! WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS CONTAINED IN THE PRESENT
! FILE. THEIR ARGUMENTS P AND Q TAKE THE VALUES L1P, L1P+/-1 AND L2P, L2P+/-1,
! RESPECTIVELY. THE BOUNDS OF THE PSIP AND PSIM ARRAYS ARE SUCH THAT ALL THE P AND Q
! VALUES REQUIRED TO PERFORM THE RECURSION ARE PROVIDED.
! THE MENTIONED "FLPHI" ROUTINE USES THE FUNCTION WIGMOD WHICH IS ALSO CONTAINED
! IN THE PRESENT FILE.
! THE PARAMETER FIRSTCALL IS A LOGICAL VARIABLE WHICH ENSURES THAT FLPHI
! IS CALLED ONLY ONCE AT A GIVEN L VALUE PROVIDING THAT ALL THE COMPUTATIONS ARE
! PERFORMED WITH THE SAME NQMAX. UNNECESSARY REPEATED CALLS OF FLPHI ARE
! SUPPRESSED VIA SAYING FIRSTCALL=.FALSE. IN A PROPER PLACE. IF THE CALLING PROGRAM
! INCLUDES NESTED LOOPS OVER N-TYPE VARIABLES FROM N=0 UP TO N=NMAX AND THESE
! LOOPS CONTAIN CALLS OF OSBRAC THEN FIRSTCALL=.FALSE. IS TO BE SAID AFTER THE MOST
! INNER OF THESE LOOPS WITH THIS PURPOSE. AND AFTER A CHANGE OF L ONE SHOULD
! CALL FLPHI ANEW VIA SAYING FIRSTCALL=.TRUE. (SEE ALSO THE APPENDED PROGRAM
! TESTOSBRAC AS AN EXAMPLE.)
!----------------------------------------------------------------------------------------------------------------------------------------
DOUBLE PRECISION BRAC,FL,PSIP,PSIM,A,FAC,DFAC,DEFAC,RFAC,BI,CO,SI,T,CO2,SI2,SC,FA,&
PLFACT,D2L12P,TN,PA,S,FA1,FA2
LOGICAL FIRSTCALL
DIMENSION BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2),&
! THIS ARRAY IS OF THE FORM: BRAC(NP,N1P,MP,N1,N2)
A(0:42,0:84),PSIP(85,85),PSIM(85,85),&
FL(0:84,0:42,0:84),&
! THIS ARRAY IS OF THE FORM: FL(NP,MP,N)
FAC(0:169),DFAC(0:84),DEFAC(0:42),RFAC(0:84),BI(0:169,0:169)
SAVE PSIP,PSIM,FL,FAC,DFAC,RFAC,A,BI
IF (NQMAX.GT.84) THEN
WRITE(6,*)'IN THE R(*8) COMPUTATION NQMAX SHOULD NOT EXCEED 84'
STOP
ENDIF
IF (CO.EQ.0.D0.OR.SI.EQ.0.D0) THEN
WRITE(6,*)'THE PROGRAM IS OF NO USE AT ZERO COS(PHI) OR SIN(PHI) VALUES'
WRITE(6,*)'COS(PHI)=',CO,' SIN(PHI)=',SI
STOP
ENDIF
IF (NQMAX.LT.L) THEN
WRITE(6,*)'NQMAX=',NQMAX,' L=',L
WRITE(6,*)'L SHOULD NOT EXCEED NQMAX'
STOP
ENDIF
! BRACKETS [N_1'L_1'N_2'L_2'|0L_10L_2]_L^\PHI TO START THE RECURSION, EQ. (22) IN
! THE ACCOMPANYING CPC PAPER.
! SUBSIDIARY QUANTITIES:
NQML=NQMAX-L
NQMLD2=NQML/2
NN=NQML-2*(NQMLD2)
NPMAX=L-NN
IF (M.GT.NQMLD2.OR.N.GT.NPMAX) THEN
WRITE(6,*) 'M SHOULD NOT EXCEED (NQ-L)2 AND N SHOULD NOT EXCEED L-NN. BUT'
WRITE(6,*) 'M=',M,' MMAX=',MPMAX,' N=',N,' NMAX=',NPMAX
STOP
ENDIF
L1=M+N+NN
L2=M-N+L
L1L2=L1+L2
L1L2ML=L1L2-L
MPMAX=(L1L2ML-NN)/2
FA=SQRT((2.D0*L1+1)*(2*L2+1))*CO**L1L2
T=SI/CO
IF (FIRSTCALL) THEN
CALL ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
CALL COE(2*NQMAX+1,FAC,BI)
! THE A_{NL} ARRAY:
DO NI=0,NQMAX/2
NA=NI-2*(NI/2)
K=1
IF(NA.EQ.1)K=-1
DO LI=0,NQMAX-2*NI
A(NI,LI)=K/SQRT(DEFAC(NI)*DFAC(NI+LI))
ENDDO
ENDDO
ENDIF
PLFACT=SQRT(FAC(L1L2+L+1)*FAC(L1L2ML))
IF (FIRSTCALL) CALL FLPHI(NQMAX,L,N,T,BI,RFAC,FL)
! THIS SUBROUTINE PRODUCES THE ARRAY FL USED BELOW.
! COMPUTATION OF THE N1=N2=0 BRACKETS. THEY ARE THE QUANTITIES BRAC(NP,N1P,MP,0,0):
DO MP=0,MPMAX
L12P=2*MP+L+NN
! L12P=L1P+L2P
D2L12P=1.D0/2.D0**L12P
N12P=MPMAX-MP
! N12P=N1P+N2P=(L1L2-L12P)/2
DO N1P=0,N12P
N1PA=N1P-2*(N1P/2)
IF (N1PA.EQ.0) THEN
K=1
ELSE
K=-1
ENDIF
! K=(-1)**N1P
TN=T**N12P*K
N2P=N12P-N1P
DO NP=0,NPMAX
L1P=MP+NP+NN
L2P=L12P-L1P
! BY DEFINITION ONE ALSO HAS: L2P=L2P(MP,NP)=MP-NP+L
PA=A(N1P,L1P)*A(N2P,L2P)
BRAC(NP,N1P,MP,0,0)=FL(NP,MP,N)*FA*TN*D2L12P*PLFACT*PA*PA
ENDDO
ENDDO
ENDDO
! RECURSION TO OBTAIN THE GENERAL FORM BRACKETS
! [N_1'L_1'N_2'L_2'|N_1L_1N_2L_2]_L^\PHI.
IF (FIRSTCALL) CALL COEFREC((NQMAX+L)/2+1,L,NN,NQMAX,PSIP,PSIM)
N12MAX=(NQMAX-L1L2)/2
MPMAX0=MPMAX
CO2=CO**2
SI2=SI**2
SC=SI*CO
! THE N1=0, N2-1-->N2 RECURSION, EQ. (20) IN THE ACCOMPANYING CPC PAPER WITH THE
! MODIFICATION POINTED OUT THERE:
DO N2=1,N12MAX
MPMAX=MPMAX+1
! THIS MPMAX VALUE IS (NQ-L-NN)/2 WHERE NQ=L1+L2+2*N2.
DO MP=0,MPMAX
N12P=MPMAX-MP
DO N1P=0,N12P
! MP=(L1P+L2P-L-NN)/2. THEREFORE, N12P=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.
DO NP=0,NPMAX
L1P=MP+NP+NN
L2P=MP-NP+L
S=0.D0
! BELOW THE RESTRICTIONS ARE IMPOSED ON THE BRAC ARRAY ENTERING THE RIGHT-HAND
! SIDE OF THE RECURRENCE FORMULAE. THIS ARRAY IS OF THE FORM BRAC(K1,K2,K3,...) WHERE
! K1=K1(NP), K2=N2(N1P), AND K3=K3(MP). NAMELY, K1=NP, OR NP+1, OR NP-1;
! K2=N1P OR N1P-1; K3=MP, OR MP-1, OR MP+1. (THE COMBINATION K2=N1P AND K3=MP+1
! DOES NOT ARISE IN THE RECURRENCE FORMULAE.) THE QUANTITIES K1, K2, AND K3
! REPRESENT, RESPECTIELY, NP, N1P, AND MP VALUES AT THE PRECEDING STAGE OF THE
! RECURSION. THE RESTRICTIONS ENSURE THAT K1, K2, AND K3 RANGE WITHIN THE LIMITS
! PERTAINING TO THE BRAC ARRAY OBTAINED AT THAT PRECEDING STAGE OF THE RECURSION.
! THUS K1, K2, AND K3 SHOULD BE NON-NEGATIVE. THEREFORE, WHEN K1=NP-1 THE VALUE
! NP=0 IS TO BE EXCLUDED, WHEN K2=N1P-1 THE VALUE N1P=0 IS TO BE EXCLUDED, AND
! WHEN K3=MP-1 THE VALUE MP=0 IS TO BE EXCLUDED.
! FURTHERMORE, IT SHOULD BE K1 \LE NMAX AND AT THE SAME TIME ONE HAS NP \LE NMAX.
! THEREFORE, WHEN K1=NP+1 THE VALUE NP=NPMAX IS TO BE EXCLUDED.
! IN ADDITION, IT SHOULD BE K2+K3 \LE MPMAX-1 WHILE ONE HAS N1P+MP \LE MPMAX. WHEN
! K2=N1P AND K3=MP-1, OR K2=N1P-1 AND K3=MP-1, OR K2=N1P-1 AND K3=MP ONE
! AUTOMATICALLY HAS N1P+MP \LE MPMAX. THEREFORE, IN THESE CASES THE CONDITION
! K2+K3 \LE MPMAX-1 DOES NOT CREATE ANY RESTRICTIONS.
! BUT WHEN K2=N1P AND K3=MP, OR K2=N1P-1 AND K3=MP+1 THE CASE N1P+MP=MPMAX
! IS TO BE FORBIDDEN.
! IN THE RECURRENCE FORMULAE BELOW THE DESCRIBED RESTRICTIONS ARE
! IMPOSED.
IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,0,N2-1)*PSIP(L1P,L2P)
IF (N1P.NE.0.AND.N1P.NE.N12P) &
S=S+BRAC(NP,N1P-1,MP+1,0,N2-1)*PSIP(L1P+1,L2P+1)
! RECALL THAT N12P=MPMAX-MP.
IF (NP.NE.NPMAX.AND.N1P.NE.0) &
S=S-BRAC(NP+1,N1P-1,MP,0,N2-1)*PSIM(L1P+1,L2P)
IF (NP.NE.0.AND.N1P.NE.N12P) &
S=S-BRAC(NP-1,N1P,MP,0,N2-1)*PSIM(L1P,L2P+1)
S=S*SC
IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,0,N2-1)*SI2
IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,0,N2-1)*CO2
BRAC(NP,N1P,MP,0,N2)=S
ENDDO
ENDDO
ENDDO
ENDDO ! N2
! N1-1-->N1 RECURSION, EQ. (20) IN THE ACCOMPANYING CPC PAPER:
DO N2=0,N12MAX
MPMAX=MPMAX0+N2
! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2=MPMAX0+N2 SINCE NQ=L1+L2+2*N2.
! THE RECURSION:
DO N1=1,N12MAX-N2
MPMAX=MPMAX+1
! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2 AND NQ=L1+L2+2*N2+2*N1.
DO MP=0,MPMAX
N12P=MPMAX-MP
DO N1P=0,N12P
! MP=(L1P+L2P-L-NN)/2. THEREFORE, N1PMAX=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.
DO NP=0,NPMAX
L1P=MP+NP+NN
L2P=MP-NP+L
S=0.D0
IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,N1-1,N2)*PSIP(L1P,L2P)
IF (N1P.NE.0.AND.N1P.NE.N12P) &
S=S+BRAC(NP,N1P-1,MP+1,N1-1,N2)*PSIP(L1P+1,L2P+1)
IF (NP.NE.NPMAX.AND.N1P.NE.0) &
S=S-BRAC(NP+1,N1P-1,MP,N1-1,N2)*PSIM(L1P+1,L2P)
IF (NP.NE.0.AND.N1P.NE.N12P) &
S=S-BRAC(NP-1,N1P,MP,N1-1,N2)*PSIM(L1P,L2P+1)
! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
S=-S*SC
IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,N1-1,N2)*CO2
IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,N1-1,N2)*SI2
! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
BRAC(NP,N1P,MP,N1,N2)=S
ENDDO
ENDDO
ENDDO
ENDDO ! N1
ENDDO ! N2
! RENORMALIZATION OF THE BRACKETS, EQ. (18) IN THE ACCOMPANYING CPC PAPER:
DO N2=0,N12MAX
MPMAX1=MPMAX0+N2
DO N1=0,N12MAX-N2
MPMAX=MPMAX1+N1
FA1=A(N1,L1)*A(N2,L2)
DO MP=0,MPMAX
N12P=MPMAX-MP
DO N1P=0,N12P
N2P=N12P-N1P
DO NP=0,NPMAX
L1P=MP+NP+NN
L2P=MP-NP+L
FA2=A(N1P,L1P)*A(N2P,L2P)
BRAC(NP,N1P,MP,N1,N2)=BRAC(NP,N1P,MP,N1,N2) &
*FA1/FA2
ENDDO ! NP
ENDDO ! N1P
ENDDO ! MP
ENDDO ! N1
ENDDO ! N2
RETURN
END
SUBROUTINE ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
! FAC(I), DFAC(I),DEFAC(I), AND RFAC(I) ARE, RESPECTIVELY, THE QUANTITIES
! I!, (2I+1)!!, (2I)!!, AND SQRT((2*I)!)/I!
DOUBLE PRECISION FAC,DFAC,DEFAC,RFAC
DIMENSION FAC(0:169),DFAC(0:84),DEFAC(0:42),RFAC(0:84)
FAC(0)=1.D0
DFAC(0)=1.D0
DEFAC(0)=1.D0
RFAC(0)=1.D0
DO I=1,2*NQMAX+1
FAC(I)=FAC(I-1)*I
ENDDO
DO I=1,NQMAX+1
DFAC(I)=DFAC(I-1)*(2*I+1)
ENDDO
DO I=1,NQMAX/2+1
DEFAC(I)=DEFAC(I-1)*2*I
ENDDO
DO I=1,NQMAX
RFAC(I)=RFAC(I-1)*2*SQRT(1-0.5D0/I)
ENDDO
RETURN
END
SUBROUTINE FLPHI(NQMAX,L,N,T,BI,RFAC,FL)
! PROVIDES THE QUANTITY SQRT((2L1P+1)*(2L2P+1))*F_L^\VARPHI, SEE EQ. (23) IN THE CPC
! PAPER, IN THE FORM OF THE FL ARRAY.
! USES THE FUNCTION WIGMOD.
! THIS FLPHI ROUTINE IS NOT IDENTICAL TO THAT IN THE OTHER FILE allosbrac.f90.
! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION BI,RFAC,FL,T,T2,SQP,F,WIGMOD,z
DIMENSION BI(0:169,0:169),RFAC(0:84),FL(0:84,0:42,0:84)
! THE OUTPUT IS FL(NP,MP,N) WHERE NP=(L1P-L2P+L-NN)/2,
! MP=(L1P+L2P-(L+NN))/2, L1P+L2P=L1T+L2T, AND N=(L1-L2+L-NN)/2, L1-L2=L1T-L2T.
T2=T*T
NQMAL=NQMAX-L
MPMAX=NQMAL/2
NN=NQMAL-2*MPMAX
LMNN=L-NN
L1ML2=2*N-LMNN
LL3=2*L
DO MP=0,MPMAX
L1PL2P=2*MP+L+NN
DO NP=0,LMNN
L1PML2P=2*NP-LMNN
L1P=(L1PL2P+L1PML2P)/2
L2P=(L1PL2P-L1PML2P)/2
M3=L2P-L1P
SQP=SQRT((2*L1P+1.D0)*(2*L2P+1))
L1T=(L1PL2P+L1ML2)/2
L2T=(L1PL2P-L1ML2)/2
N3=L1T+L2T-L
IMIN=ABS(L1T-L1P)
IMAX=MIN(L1P+L1T,L2P+L2T)
NAL1=(L1T+L1P-IMIN)/2
NAL2=(L1T-L1P+IMIN)/2
NAL3=(L2T-L2P+IMIN)/2
NAL4=(L2T+L2P-IMIN)/2
NAL4P=NAL4-2*(NAL4/2)
F=0.D0
z=t**imin
if(nal4p.ne.0)z=-z
DO I=IMIN,IMAX,2
F=F+RFAC(NAL1)*RFAC(NAL2)*RFAC(NAL3)*RFAC(NAL4)*&
WIGMOD(L1T,L2T,L,L1P-I,I-L2P,BI)*z
NAL1=NAL1-1
NAL2=NAL2+1
NAL3=NAL3+1
NAL4=NAL4-1
z=-z*t2
ENDDO
FL(NP,MP,N)=F*SQRT((2*L1P+1)*(2*L2P+1)&
*BI(L1T+L-L2T,2*L1T)*BI(N3,2*L2T)/((LL3+1)*BI(N3,L1T+L2T+L+1)*BI(L+M3,LL3)))
ENDDO
ENDDO
RETURN
END
FUNCTION WIGMOD(L1,L2,L3,M1,M2,BI)
! THE 3J SYMBOL IN TERMS OF BINOMIAL COEFFICIENTS WITHOUT THE FACTOR
! SQRT(F) WHERE F IS AS FOLLOWS,
! F=BI(L1+L3-L2,2*L1)*BI(N3,2*L2)/((LL3+1)*BI(N3,L1+L2+L3+1)*BI(L3+M3,LL3))
! WITH N3=L1+L2-L3, LL3=2*L3, AND M3=-M1-M2.
! THIS WIGMOD ROUTINE IS NOT IDENTICAL TO THAT IN THE OTHER FILE allosbrac.f90.
DOUBLE PRECISION WIGMOD,BI,S
DIMENSION BI(0:169,0:169)
M3=-M1-M2
N3=L1+L2-L3
LM1=L1-M1
LP2=L2+M2
KMIN=MAX(0,L1-L3+M2,L2-L3-M1)
KMAX=MIN(N3,LM1,LP2)
S=0.d0
NPH=1
DO K=KMIN,KMAX
S=S+NPH*BI(K,N3)*BI(LM1-K,L1+L3-L2)*BI(LP2-K,L2+L3-L1)
NPH=-NPH
ENDDO
WIGMOD=S/SQRT(BI(LM1,2*L1)*BI(LP2,2*L2))
NY=KMIN+L1-L2-M3
NYP=NY-2*(NY/2)
IF(NYP.NE.0)WIGMOD=-WIGMOD
RETURN
END
SUBROUTINE COE(NMAX,FAC,BI)
! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION FAC,BI
DIMENSION FAC(0:169),BI(0:169,0:169)
DO N=0,NMAX
DO M=0,N/2
BI(M,N)=FAC(N)/(FAC(M)*FAC(N-M))
BI(N-M,N)=BI(M,N)
ENDDO
ENDDO
RETURN
END
SUBROUTINE COEFREC(MAXCOE,L,NN,NQMAX,PSIP,PSIM)
! CALCULATES THE COEFFICIENTS OF THE RECURSION FORMULA, EQ. (21) IN THE ACCOMPANYING
! CPC PAPER
DOUBLE PRECISION PSIP,PSIM
DIMENSION PSIP(85,85),PSIM(85,85)
! MAXCOE EQUALS INT((NQMAX+LMAX)/2)+1 AT CALLS OF THIS ROUTINE.
LP=L+1
LM=L-1
DO M1=1,MAXCOE
DO M2=1,M1
M1P2=M1+M2
NA=M1P2+L
NNA=NA-2*(NA/2)
M1M2=M1-M2
! CALCULATION OF PSIP. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P AND L2P OR
! L1P+1 AND L2P+1.)
IF (NNA.EQ.NN.AND.M1P2.GE.L.AND.M1P2.LE.NQMAX.AND.&
ABS(M1M2).LE.L) THEN
IF (M1P2.GT.LP) THEN
M1P2L=M1P2+L
M1P2ML=M1P2-L
PSIP(M1,M2)=SQRT(M1P2L*(M1P2L+1.D0)*M1P2ML*(M1P2ML-1)/((4*M1*M1-1)&
*(4*M2*M2-1)))
ELSE
PSIP(M1,M2)=0.D0
ENDIF
PSIP(M2,M1)=PSIP(M1,M2)
ENDIF
! CALCULATION OF PSIM. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P+1 AND L2P
! OR L2P+1 AND L1P.)
IF (NNA.NE.NN.AND.M1P2.GE.LP.AND.M1P2.LE.NQMAX&
.AND.M1M2.LE.LP.AND.M1M2.GE.-LM) THEN
IF (M1M2.LT.L) THEN
M1M2L=M1M2+L
M1M2ML=M1M2-L
PSIM(M1,M2)=SQRT(M1M2L*(M1M2L+1.D0)*M1M2ML*(M1M2ML-1)/&
((4*M1*M1-1)*(4*M2*M2-1)))
ELSE
PSIM(M1,M2)=0.D0
ENDIF
PSIM(M2,M1)=PSIM(M1,M2)
ENDIF
ENDDO
ENDDO
RETURN
END

260
OSBRACKETS/testallosbrac.f90
View File

@ -1,260 +0,0 @@
program testallosbrac
! Tests of the program that calculates the <N_1'L_1'N_2'L_2'|N_1L_1N_2L_2>_L^\varphi brackets.
implicit double precision(a-h,o-z)
dimension fac(0:201),dfac(0:201),defac(0:201),rfac(0:100)
! These arrays enter the routine ARR called below. The first three arrays are used merely to
! compute the exact expression for the bracket below. The array rfac is of no use.
allocatable :: brac(:,:,:,:,:,:,:,:),bi(:,:)
! Below the allocation/deallocation actions are done independently for each of the tests
! performed as it would be required in the case when only a single test is being performed.
open(13,file='allosoutput')
! The input angular parameters are CO=COS(VARPHI) and SI=SIN(VARPHI). Set, for example,
si=1/sqrt(2.d0)
co=-si
! This corresponds to Eq. (13) in the accompanying CPC paper in the equal mass case. The
! corresponding transformation (7) in that paper arises at calculating standard two-body
! shell-model matrix elements.
! One might also test the program for COS(VARPHI) and SIN(VARPHI) values pertaining to
! nonequal mass cases. For example, one might set
! si=1/sqrt(3.d0)
! co=-sqrt(2/3.d0)
! 1. Comparison of the value of the <n1p,L00|n1,l1,n2,l2>_L^\varphi bracket with the result
! of its analytic calculation. Here l1+l2 should be of the same parity as L.
l=4
l1=2
l2=4
n1=0
n2=0
nq=2*n1+2*n2+l1+l2
! also nq=2*n1p+L
NN=0
! this is because nq-l=2*n1p is even
mp=0
! here l1p=l, mp=(l1p+l2p-l-nn)/2=0
np=l
! np=(l1p-l2p+l-nn)/2=l
n1p=(nq-l)/2
! Computation of the exact bracket.
nfac=2*nq+1
CALL ARR(FAC,DFAC,DEFAC,RFAC,NQ)
! The array of binomial coefficients to calculate the Wigner coefficients:
allocate (BI(0:NFAC,0:NFAC),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate BI failed"
call coe(NFAC,FAC,BI)
nm=n1p-n1-n2
nfa=nm-2*(nm/2)
LL3=2*L
N3=L1+L2-L
! exact bracket:
b=co**(2*n1+l1)*si**(2*n2+l2)*(-1)**l*sqrt((2.d0*l1+1)*(2*l2+1))*(-1)**nfa*&
WIGMOD(l1,l2,l,0,0,bi,nfac)*SQRT(BI(L1+L-L2,2*L1)*BI(N3,2*L2)/&
((LL3+1)*BI(N3,L1+L2+L+1)*BI(L,LL3)))&
*sqrt(defac(n1p)*dfac(n1p+l)/(defac(n1)*defac(n2)*dfac(n1+l1)*dfac(n2+l2)))
! bracket from the code:
ndel=4
nqmax=nq+ndel
! The result should not depend on ndel, provided that ndel is even, i.e., nq and nqmax are of the
! same parity.
NQML=NQMAX-L
m=(l1+l2-l-nn)/2
n=(l1-l2+l-nn)/2
allocate (BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,&
0:(NQMAX-L)/2,0:L,0:(NQMAX-L)/2,L:L),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
CALL ALLOSBRAC(NQMAX,L,L,CO,SI,BRAC)
a=brac(np,n1p,mp,n1,n2,n,m,l)
write(6,*)'1. bracket <n1p L 0 0|n1 l1 n2 l2>_L^varphi'
write(6,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l
write(6,*)'exact bracket',b,' bracket from the code',a
write(13,*)'1. bracket <n1p L 0 0|n1 l1 n2 l2>_L^varphi'
write(13,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l
write(13,*)'exact bracket',b,' bracket from the code',a
! write(6,*)co,si
deallocate (BRAC,STAT=istatus)
! 2. Test of the symmetry of the brackets. The relation
! <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\varphi=<N1,L1,N2,L2|N1P,L1P,N2P,L2P>_L^\varphi
! should hold true. This symmetry relation may be rewritten as
! BRAC(NP,N1P,MP,N1,N2,N,M,L)=BRAC1(N,N1,M,N1P,N2P,NP,MP,L) where N and
! M represent the angular momenta L1 and L2 and NP and MP represent L1P and L2P.
! An example:
nq=30
l=9
! nn=1 in the present example.
! calculate: nmax=l-nn.
! nmax=8 in the present example. Any n and np that do not exceed nmax may be employed.
n=4
np=3
! calculate: nqmld2=(nq-l-nn)/2.
! nqmld2=10 in the present example. Any m, n1, n2, mp, n1p, and n2p may be employed
! provided that m+n1+n2=mp+n1p+n2p=nqmld2.
m=5
n1=3
n2=2
mp=4
n1p=3
n2p=3
nqmax=32
! Recall that nq=30. Any nqmax of the same parity as nq and such that nqmax.ge.nq
! may be employed. Then its choice does not influence the results.
allocate (BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,&
0:(NQMAX-L)/2,0:L,0:(NQMAX-L)/2,L:L),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
CALL ALLOSBRAC(NQMAX,L,L,CO,SI,BRAC)
write(6,*)'2. Symmetry test, L=',l
write(6,*)'m =',m,' n =',n,' n1 =',n1,' n2 =',n2
write(6,*)'mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
write(6,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =',BRAC(np,n1p,mp,n1,n2,n,m,l)
write(6,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =',BRAC(n,n1,m,n1p,n2p,np,mp,l)
write(13,*)'2. Symmetry test, L=',l
write(13,*)'m =',m,' n =',n,' n1 =',n1,' n2 =',n2
write(13,*)'mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
write(13,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =',BRAC(np,n1p,mp,n1,n2,n,m,l)
write(13,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =',BRAC(n,n1,m,n1p,n2p,np,mp,l)
deallocate (BRAC,STAT=istatus)
! 3. Test of the relation
! \sum_{i,i',L}\sum_j<j|i><j|i'>=\sum_L N_0(L).
! Here i=(n_1,l_1,n_2,l_2), i'=(n_3,l_3,n_4,l_4), and j=(n_1p,l_1p,n_2p,l_2p). The l.h.s. sum
! runs over all i and i' values such that nq(i)=nq(i')=nq.
! In the r.h.s. sum, N_0(L) is the number of the (n_1,l_1,n_2,l_2) states pertaining to
! given nq and L values. The sum runs over all the L values compatible with a given nq value
! and thus represents the total number of states with this nq.
nq=14
ndel=4
nqmax=nq+ndel
! The result should not depend on ndel, i.e., on nqmax, provided that ndel is even.
numb=0
s=0.d0
LMIN=0
LMAX=NQ
allocate (BRAC(0:LMAX,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,&
0:(NQMAX-LMIN)/2,0:LMAX,0:(NQMAX-LMIN)/2,LMIN:LMAX),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
CALL ALLOSBRAC(NQMAX,LMIN,LMAX,CO,SI,BRAC)
do L=0,nq
nqml=nq-l
nqmld2=nqml/2
nn=nqml-2*nqmld2
mmax=(nq-(l+nn))/2
nmax=l-nn
nst=(nmax+1)*(mmax+1)*(mmax+2)/2
do na=0,nmax
do ma=0,mmax
! IN THE ABOVE NOTATION NA=NA(L1,L2), MA=MA(L1,L2).
n12=mmax-ma
do nb=0,nmax
do mb=0,mmax
! IN THE ABOVE NOTATION NB=NB(L3,L4), MB=MB(L4,L4).
n34=mmax-mb
ds=0.d0
do n2=0,n12
n1=n12-n2
do n4=0,n34
n3=n34-n4
DO MP=0,mmax
do n1p=0,mmax-mp
DO NP=0,nmax
DS=DS+BRAC(np,n1p,mp,n1,n2,na,ma,l)&
*BRAC(np,n1p,mp,n3,n4,nb,mb,l)
enddo ! NP
ENDDO ! n1p
ENDDO ! MP
enddo ! n4
enddo ! n2
s=s+ds
enddo ! MB
enddo ! NB
enddo ! MA
enddo ! NA
numb=numb+nst
ENDDO !l
deallocate (BRAC,STAT=istatus)
write(6,*)'3. Test of the relation sum_{i,ip,L}sum_j<j|i><j|ip> = the number of states.'
write(6,*)'See the comments in the text of the program and in more detail in the'
write(6,*)'accompanying CPC paper for further explanations.'
write(6,*)'nqmax =',nqmax,' nq =',nq
write(6,*)'sum_{i1,i2,L}sum_j<j|i1><j|i2> =',s,' exact value (number of states) =', numb
write(13,*)'3. Test of the relation sum_{i,ip,L}sum_j<j|i><j|ip> = the number of states.'
write(13,*)'See the comments in the text of the program and in more detail in the'
write(13,*)'accompanying CPC paper for further explanations.'
write(13,*)'nqmax =',nqmax,' nq =',nq
write(13,*)'sum_{i1,i2,L}sum_j<j|i1><j|i2> =',s,' exact value (number of states)=',numb
! 4. Test of the relation \sum_{i,i',L}|\sum_j<j|i><j|i'>-\delta_{i,i'}|=0.
! Here i=(n_1,l_1,n_2,l_2), i'=(n_1',l_1',n_2',l_2'), and j=(n_1'',l_1'',n_2'',l_2''). The outer
! sum runs over all i and i' values such that nq(i)=nq(i') \le nqmax.
! Also one has nq(i)=nq(i')=nq(j).
nqmax=12
lmin=0
lmax=nqmax
nqml=nqmax-lmin
nqmld2=nqml/2
allocate (BRAC(0:LMAX,0:nqmld2,0:nqmld2,0:nqmld2,0:nqmld2,0:LMAX,&
0:nqmld2,LMIN:LMAX),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
CALL ALLOSBRAC(NQMAX,LMIN,LMAX,CO,SI,BRAC)
! When nqmax is sufficiently high, this calculation is to be performed in somewhat another way
! because of the memory restrictions. In such cases one may use ALLOSBRAC at lmin=lmax=l.
! The BRAC array is to be allocated and deallocated inside the loop over l below, and
! ALLOSBRAC is to be called inside this loop.
s=0.d0
! The quantity s represents the expression \sum_{i,i',L}|\sum_j<j|i><j|i'>-\delta_{i,i'}| we are
! calculating.
do l=0,nqmax
nqml=nqmax-l
nqmld2=nqml/2
nn=nqml-2*nqmld2
nmax=l-nn
lpnn=l+nn
do nq=lpnn,nqmax,2
mmax=(nq-lpnn)/2
do ma=0,mmax
! The equality mmax=(nq-lpnn)/2 follows from the definition of m-type variables.
n12=mmax-ma
! The quantity n12 designates n1+n2 and here it is taken into account that
! n1+n2+ma=mmax.
do na=0,nmax
! In the case of odd nqmax and l equal to zero, one has nmax=-1 and no contribution arises as it
! should be.
do n1=0,n12
n2=n12-n1
do mb=0,mmax
n34=mmax-mb
do nb=0,nmax
do n3=0,n34
n4=n34-n3
ds=0.d0
! The quantity ds represents the contribution of \sum_j|<j|i>_L<j|i'>_L-\delta_{i,i'}| to the net
! result.
DO MP=0,MMAX
do n1p=0,mmax-mp
DO NP=0,NMAX
DS=DS+BRAC(np,n1p,mp,n1,n2,na,ma,l)&
*BRAC(np,n1p,mp,n3,n4,nb,mb,l)
enddo ! NP
ENDDO ! n1p
ENDDO ! MP
if(ma.eq.mb.and.na.eq.nb.and.n1.eq.n3.and.n2.eq.n4)then
ds=abs(ds-1.d0)
else
ds=abs(ds)
endif
s=s+ds
enddo ! n3
enddo ! nb
enddo ! mb
enddo ! n1
enddo ! na
enddo ! ma
enddo ! nq
ENDDO ! l
deallocate (BRAC,STAT=istatus)
write(6,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
write(6,*)'For further explanations see the comments in the text of the program'
write(6,*)'and, in more detail, in the accompanying CPC paper.'
write(6,*)'nqmax =',nqmax
write(6,*)'sum_{i1,i2,L}|sum_j<j|i1><j|i2>-delta_{i1,i2}| =',s
write(13,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
write(13,*)'For further explanations see the comments in the text of the program'
write(13,*)'and, in more detail, in the accompanying CPC paper.'
write(13,*)'nqmax =',nqmax
write(13,*)'sum_{i1,i2,L}|sum_j<j|i1><j|i2>-delta_{i1,i2}| =',s
end

266
OSBRACKETS/testosbrac.f90
View File

@ -1,266 +0,0 @@
program testosbrac
! Tests of the program that calculates the <N_1'L_1'N_2'L_2'|N_1L_1N_2L_2>_L^\varphi brackets.
implicit double precision(a-h,o-z)
logical firstcall
dimension fac(0:169),dfac(0:84),defac(0:42),rfac(0:84),bi(0:169,0:169)
! These arrays enter the routine ARR called below. The first three arrays are used merely to
! compute the exact expression for the bracket below. The array rfac is of no use.
allocatable :: brac(:,:,:,:,:),brac1(:,:,:,:,:)
open(13,file='osoutput')
! The input angular parameters are CO=COS(VARPHI) and SI=SIN(VARPHI). Set, for example,
si=1/sqrt(2.d0)
co=-si
! This corresponds to Eq. (13) in the accompanying CPC paper in the equal mass case. The
! corresponding transformation (7) in that paper arises at calculating standard two-body
! shell-model matrix elements.
! One might also test the program for COS(VARPHI) and SIN(VARPHI) values pertaining to
! nonequal mass cases. For example, one might set
! si=1/sqrt(3.d0)
! co=-sqrt(2/3.d0)
! 1. Comparison of the value of the <n1p,L00|n1,l1,n2,l2>_L^\varphi bracket with the result
! of its analytic calculation. Here l1+l2 should be of the same parity as L.
l=4
l1=2
l2=4
n1=3
n2=2
mp=0
np=l
nq=2*n1+2*n2+l1+l2
n1p=(nq-l)/2
! Computation of the exact bracket.
nfac=2*nq+1
CALL ARR(FAC,DFAC,DEFAC,RFAC,NQ)
! The array of binomial coefficients to calculate the Wigner coefficients:
call coe(NFAC,FAC,BI)
nm=n1p-n1-n2
nfa=nm-2*(nm/2)
LL3=2*L
N3=L1+L2-L
! exact bracket:
b=co**(2*n1+l1)*si**(2*n2+l2)*(-1)**l*sqrt((2.d0*l1+1)*(2*l2+1))*(-1)**nfa*&
WIGMOD(l1,l2,l,0,0,bi)*SQRT(BI(L1+L-L2,2*L1)*BI(N3,2*L2)/&
((LL3+1)*BI(N3,L1+L2+L+1)*BI(L,LL3)))&
*sqrt(defac(n1p)*dfac(n1p+l)/(defac(n1)*defac(n2)*dfac(n1+l1)*dfac(n2+l2)))
! bracket from the code:
ndel=4
nqmax=nq+ndel
! The result should not depend on ndel, provided that ndel is even, i.e., nq and nqmax are of the
! same parity.
NQML=NQMAX-L
NN=0
m=(l1+l2-l-nn)/2
n=(l1-l2+l-nn)/2
allocate (BRAC(0:L,0:NQML/2,0:NQML/2,0:NQML/2,0:NQML/2),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
FIRSTCALL=.TRUE.
CALL OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
a=brac(np,n1p,mp,n1,n2)
write(6,*)'1. bracket <n1p,L00|n1l1n2l2>_L^phi'
write(6,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l
write(13,*)'1. bracket <n1p,L00|n1l1n2l2>_L^phi'
write(13,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l
write(6,*)'exact bracket',b,' bracket from the code',a
write(13,*)'exact bracket',b,' bracket from the code',a
deallocate (BRAC,STAT=istatus)
! 2. Test of the symmetry of the brackets. The relation
! <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\varphi=<N1,L1,N2,L2|N1P,L1P,N2P,L2P>_L^\varphi
! should hold true. Let the bracket in the left-hand side of this relation is contained in the BRAC
! array and that in its right-hand side in the BRAC1 array. Then the above symmetry relation may
! be rewritten as BRAC(NP,N1P,MP,N1,N2)=BRAC1(N,N1,M,N1P,N2P) where N and M represent
! the angular momenta L1 and L2 and NP and MP represent L1P and L2P.
! An example:
nq=30
l=9
! nn=1 at this choice.
! calculate: nmax=l-nn.
! nmax=8 in the present example. Any n and np values
! that do not exceed nmax may be employed.
n=4
np=3
! calculate: nqmld2=(nq-l-nn)/2.
! nqmld2=10 in the present example. Any m, n1, n2, mp, n1p, and n2p may be employed
! provided that m+n1+n2=mp+n1p+n2p=nqmld2.
m=5
n1=3
n2=2
mp=4
n1p=3
n2p=3
nqmax=36
! Recall that nq=30. Any nqmax of the same parity as nq and such that nqmax.ge.nq
! may be employed. Then its choice does not influence the results.
nqml=nqmax-l
allocate (BRAC(0:L,0:NQML/2,0:NQML/2,0:NQML/2,0:NQML/2),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
CALL OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
allocate (BRAC1(0:L,0:NQML/2,0:NQML/2,0:NQML/2,0:NQML/2),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC1) failed"
CALL OSBRAC(NP,MP,L,NQMAX,CO,SI,FIRSTCALL,BRAC1)
write(6,*)'2. Symmetry test, L =',l
write(6,*)'m =',m,' n =',n,' n1=',n1,' n2=',n2,' mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
write(6,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =', BRAC(np,n1p,mp,n1,n2)
write(6,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =', BRAC1(n,n1,m,n1p,n2p)
write(13,*)'2. Symmetry test, L=',l
write(13,*)'m =',m,' n =',n,' n1=',n1,' n2=',n2,' mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
write(13,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =', BRAC(np,n1p,mp,n1,n2)
write(13,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =', BRAC1(n,n1,m,n1p,n2p)
deallocate (BRAC1,STAT=istatus)
deallocate (BRAC,STAT=istatus)
! 3. Test of the relation
! \sum_{i,i',L}\sum_j<j|i>_L<j|i'>_L=\sum_L N_0(L).
! Here i=(n_1,l_1,n_2,l_2), i'=(n_3,l_3,n_4,l_4), and j=(n_1p,l_1p,n_2p,l_2p). The l.h.s. sum
! runs over all i and i' values such that nq(i)=nq(i')=nq.
! In the r.h.s. sum, N_0(L) is the number of the (n_1,l_1,n_2,l_2) states pertaining to
! given nq and L values. The sum runs over all the L values compatible with a given nq value
! and thus represents the total number of states with this nq.
nq=14
ndel=6
nqmax=nq+ndel
! The result should not depend on ndel, i.e., on nqmax, provided that ndel is even.
numb=0
s=0.d0
do l=0,nq
FIRSTCALL=.true.
nqml=nqmax-l
nqmld2=nqml/2
nn=nqml-2*nqmld2
mmax=(nq-(l+nn))/2
nmax=l-nn
nst=(nmax+1)*(mmax+1)*(mmax+2)/2
allocate (BRAC(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
allocate (BRAC1(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC1) failed"
! The allocation/deallocation of BRAC and BRAC1, at nqmld2-->nqmax/2, might also be done
! outside the loop.
do ma=0,mmax
n12=mmax-ma
do na=0,nmax
! IN THE ABOVE NOTATION NA=NA(L1,L2) AND MA=MA(L1,L2).
CALL OSBRAC(NA,MA,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
! The order of summations below aims to avoid unnecessary calls of osbrac.
do mb=0,mmax
n34=mmax-mb
do nb=0,nmax
! IN THE ABOVE NOTATION NB=NB(L3,L4) AND MB=MB(L3,L4).
CALL OSBRAC(NB,MB,L,NQMAX,CO,SI,FIRSTCALL,BRAC1)
ds=0.d0
do n2=0,n12
n1=n12-n2
do n4=0,n34
n3=n34-n4
DO MP=0,mmax
do n1p=0,mmax-mp
DO NP=0,nmax
DS=DS+BRAC(NP,n1p,MP,n1,n2)*BRAC1(NP,n1p,MP,n3,n4)
enddo ! NP
ENDDO ! n1p
ENDDO ! MP
enddo ! n4
enddo ! n2
s=s+ds
enddo ! nb
FIRSTCALL=.false.
! If the nb loop is nested inside the mb loop the FIRSTCALL=.false. command is more efficient.
! In practice, usually the gain in the running time is relatively not large.
enddo ! mb
enddo ! ma
enddo ! na
numb=numb+nst
deallocate (BRAC,STAT=istatus)
deallocate (BRAC1,STAT=istatus)
ENDDO ! l
write(6,*)'3. Test of the relation sum_{i,ip,L,j}<j|i>_L<j|ip>_L = the number of states.'
write(6,*)'See the comments in the text of the program and in more detail in the'
write(6,*)'accompanying CPC paper for further explanations.'
write(6,*)'nqmax =',nqmax,' nq =',nq
write(6,*)'sum_{i,ip,L,j}<j|i>_L<j|ip>_L =',s,' exact value (number of states) =',numb
write(13,*)'3. Test of the relation sum_{i,ip,L}sum_j<j|i><j|ip> = the number of states.'
write(13,*)'See the comments in the text of the program and in more detail in the'
write(13,*)'accompanying CPC paper for further explanations.'
write(13,*)'nqmax =',nqmax,' nq =',nq
write(13,*)'sum_{i,ip,L,j}<j|i>_L<j|ip>_L =',s,' exact value (number of states) =',numb
! 4. Test of the relation \sum_{i,i',L}\sum_j|<j|i>_L<j|i'>_L-\delta_{i,i'}|=0.
! Here i=(n_1,l_1,n_2,l_2), i'=(n_1',l_1',n_2',l_2'), and j=(n_1'',l_1'',n_2'',l_2''). The outer
! sum runs over all i and i' values such that nq(i)=nq(i') \le nqmax.
! Also one has nq(i)=nq(i')=nq(j).
nqmax=12
s=0.d0
! The quantity s represents the expression \sum_{i,i',L}|\sum_j<j|i><j|i'>-\delta_{i,i'}| we are
! calculating.
do l=0,nqmax
FIRSTCALL=.true.
nqml=nqmax-l
nqmld2=nqml/2
nn=nqml-2*nqmld2
nmax=l-nn
lpnn=l+nn
allocate (BRAC(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
allocate (BRAC1(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC1) failed"
! The allocation/deallocation of BRAC and BRAC1, at nqmld2-->nqmax/2, might also be done
! outside the loop.
do ma=0,nqmld2
! The above upper limit of summation corresponds to the definition of m-type
! variables.
do na=0,nmax
! In the case of odd nqmax and l equal to zero, one has nmax=-1 and no contribution arises as it
! should be.
CALL OSBRAC(NA,MA,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
! The order of summations below aims to avoid unnecessary calls of osbrac.
do mb=0,nqmld2
do nb=0,nmax
CALL OSBRAC(NB,MB,L,NQMAX,CO,SI,FIRSTCALL,BRAC1)
do nq=l+nn,nqmax,2
nqcmld2=(nq-l-nn)/2
n12=nqcmld2-ma
n34=nqcmld2-mb
! We consider nq to be the number of quanta pertaining to states entering the brackets in our
! sum. Then one has: ma+n1+n2=mb+n3+n4=mp+n1p+n2p=(nq-l-nn)/2 that is
! nqcmld2.
do n1=0,n12
n2=n12-n1
do n3=0,n34
n4=n34-n3
ds=0.d0
! The quantity ds represents the contribution of \sum_j|<j|i>_L<j|i'>_L-\delta_{i,i'}| to the net
! result.
DO MP=0,nqcmld2
do n1p=0,nqcmld2-mp
DO NP=0,NMAX
DS=DS+BRAC(NP,n1p,MP,n1,n2)*BRAC1(NP,n1p,MP,n3,n4)
enddo ! NP
ENDDO ! n1p
ENDDO ! MP
if(ma.eq.mb.and.na.eq.nb.and.n1.eq.n3.and.n2.eq.n4)then
ds=abs(ds-1.d0)
else
ds=abs(ds)
endif
s=s+ds
enddo ! n3
enddo ! n1
enddo ! nq
enddo ! nb
FIRSTCALL=.false.
! If the nb loop is nested inside the mb loop the FIRSTCALL=.false. command is more efficient.
! In practice, usually the gain in the running time is relatively not large.
enddo ! mb
enddo ! ma
enddo ! na
deallocate (BRAC,STAT=istatus)
deallocate (BRAC1,STAT=istatus)
ENDDO ! l
write(6,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
write(6,*)'For further explanations see the comments in the text of the program'
write(6,*)'and, in more detail, in the accompanying CPC paper.'
write(6,*)'nqmax =',nqmax
write(6,*)'sum_{i1,i2,L}|sum_j<j|i1>_L<j|i2>_L-delta_{i1,i2}| =',s
write(13,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
write(13,*)'For further explanations see the comments in the text of the program'
write(13,*)'and, in more detail, in the accompanying CPC paper.'
write(13,*)'nqmax =',nqmax
write(13,*)'sum_{i1,i2,L}|sum_j<j|i1>_L<j|i2>_L-delta_{i1,i2}| =',s
end