SUBROUTINE sfplot_ABM_main(Nbasis, Nbasiss, Neigmin, rho, nang, ithrho, MegaVal) USE Numeric_Kinds_Module USE Numbers_Module USE FileUnits_Module USE Narran_Module USE Masses_Module USE chiang_Module USE Parms_Module USE convrsns_Module USE converge_Module USE totj_Module USE quantb_Module USE gaussb_Module USE InputFile_Module IMPLICIT NONE CHARACTER(LEN=5) RhoName, MegaName CHARACTER(LEN=5), PARAMETER:: blank=' ' INTEGER Plot_sf INTEGER itheta, ichi, ienergy, NTheta_Plot, nang, maxn, maxl, NChi_Plot, scheme, nthetaChi_Plot, nsfn, iang, i, j, iarran, vib INTEGER :: nbasis(0:mxmega), nbasiss(0:mxmega), neigmin(0:mxmega), ithmode, ithrho, MegaVal, NABM_Sfunc, NSYM_States, NFull_States LOGICAL :: primitive, Plot_Sfunc LOGICAL :: there REAL(KIND=SP_Kind) SP_PSI ! Single Precision for Plotting REAL(KIND=WP_Kind) u, v, dchi, dtheta, energy(10) REAL(KIND=WP_Kind) rho, theta, chi, phi, valmax, a, b, c, phi_s, rhoabm, SignATMax REAL(KIND=WP_Kind), ALLOCATABLE:: thetaval(:), chivals(:), Phi_Plot(:), Phi_ABM(:,:) REAL(KIND=WP_Kind), ALLOCATABLE:: vect(:,:), thetaonly(:), chionly(:) REAL(KIND=WP_Kind), ALLOCATABLE:: Theta_Plot(:), Chi_Plot(:), weight(:) REAL(KIND=WP_Kind), ALLOCATABLE:: phinow(:,:), sfun(:,:) REAL(KIND=WP_Kind), ALLOCATABLE:: pn(:,:,:), dpndx(:,:,:) REAL(KIND=WP_Kind), ALLOCATABLE:: oscil(:,:,:) REAL(KIND=WP_Kind), ALLOCATABLE:: doscdz(:,:,:), parityfn(:,:) REAL(KIND=WP_Kind), ALLOCATABLE:: bfaci(:,:), eigen(:,:) DATA energy/0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00/ NAMELIST/sfplot_ABM/ NTheta_Plot, a, b, c, phi_s, ithmode, primitive, Plot_Sfunc rhoabm=rho INQUIRE(file=InputDIR(1:LEN(TRIM(InputDIR)))//InputFile,exist=there) IF(there)THEN OPEN(Unit=In_Unit,File=InputDIR(1:LEN(TRIM(InputDIR)))//InputFile,status='old') READ(In_Unit,NML=sfplot_ABM, IOSTAT=plot_sf) IF(plot_sf < 0) Plot_SFunc=.False. NChi_Plot = 8*NTheta_Plot - 7 nthetaChi_Plot=NTheta_Plot*NChi_Plot nsfn=NBasis(0) rho=rhoabm CLOSE(In_Unit) IF(.NOT.Plot_Sfunc)RETURN ELSE WRITE(Msg_Unit,*)InputDIR(1:LEN(TRIM(InputDIR)))//InputFile, " does not exist." Plot_SFunc=.False. RETURN ENDIF NABM_Sfunc=NEigMin(MegaVal) NSYM_States=NBasiss(MegaVal) NFull_States=NBasis(MegaVal) IF(.NOT.ALLOCATED(ThetaOnly))ALLOCATE(ThetaOnly(NTheta_Plot)) IF(.NOT.ALLOCATED(ChiOnly))ALLOCATE(ChiOnly(NChi_Plot)) IF(.NOT.ALLOCATED(ThetaVal))ALLOCATE(thetaval(nang)) IF(.NOT.ALLOCATED(ChiVals))ALLOCATE(chivals(nang)) IF(.NOT.ALLOCATED(Weight))ALLOCATE(weight(nang)) IF(.NOT.ALLOCATED(Phi_Plot))ALLOCATE(Phi_Plot(nthetaChi_Plot)) IF(.NOT.ALLOCATED(Phi_ABM))ALLOCATE(Phi_ABM(nang,NFull_States)) IF(.NOT.ALLOCATED(PhiNow))ALLOCATE(phinow(nthetaChi_Plot,NFull_States)) IF(.NOT.ALLOCATED(SFUN))ALLOCATE(sfun(nthetaChi_Plot,NABM_Sfunc)) IF(.NOT.ALLOCATED(ParityFn))ALLOCATE(parityfn(nthetaChi_Plot,narran)) IF(.NOT.ALLOCATED(Bfaci))ALLOCATE(bfaci(nthetaChi_Plot,narran)) IF(.NOT.ALLOCATED(pn))ALLOCATE(pn(0:mxl,nthetaChi_Plot,narran)) IF(.NOT.ALLOCATED(dpndx))ALLOCATE(dpndx(0:mxl,nthetaChi_Plot,narran)) IF(.NOT.ALLOCATED(oscil))ALLOCATE(oscil(0:mxl,nthetaChi_Plot,narran)) IF(.NOT.ALLOCATED(doscdz))ALLOCATE(doscdz(0:mxl,nthetaChi_Plot,narran)) IF(.NOT.ALLOCATED(eigen))ALLOCATE(eigen(nsfn,0:mxmega)) IF(.NOT.ALLOCATED(vect))ALLOCATE(vect(NSym_States,NABM_Sfunc)) IF(.NOT.ALLOCATED(Theta_Plot))ALLOCATE(Theta_Plot(nthetaChi_Plot)) IF(.NOT.ALLOCATED(Chi_Plot))ALLOCATE(Chi_Plot(nthetaChi_Plot)) rhoabm = rho INQUIRE(FILE='plot.vecs',EXIST=there) IF(.NOT.there)THEN OPEN(Unit=PlotVecs_Unit,File=OutDIR(1:LEN(TRIM(OutDIR)))//'SfunctionPlots/PlotVecs.txt') CLOSE(UNIT=PlotVecs_Unit) ENDIF ! ! OPEN output files ! WRITE(Out_Unit,*)'Called SFPlot_ABM_Main' OPEN(Unit=Out_Unit,File=OutDIR(1:LEN(TRIM(OutDIR)))//'Output/SFPlot_ABM_Main.txt') OPEN(Unit=Plot_txt_Unit,File=OutDIR(1:LEN(TRIM(OutDIR)))//'SfunctionPlots/plot.txt') OPEN(Unit=MathPoly_Unit,File=OutDIR(1:LEN(TRIM(OutDIR)))//'SfunctionPlots/Math.Poly') OPEN(Unit=Spin_Unit,File=OutDIR(1:LEN(TRIM(OutDIR)))//'SfunctionPlots/Spin.txt') rho=rhoabm IF(abs(rho-rhoabm)>1.E-7)THEN WRITE(Out_Unit,*)'Error rhoabm/=rho' WRITE(Out_Unit,*)'rhoabm,rho', rhoabm, rho WRITE(Msg_Unit,*)'Error rhoabm/=rho' WRITE(Msg_Unit,*)'rhoabm,rho', rhoabm, rho STOP "SFPlot_ABM" ENDIF OPEN(Unit=In_Unit,File=InputDIR(1:LEN(TRIM(InputDIR)))//InputFile,status='old') READ(In_Unit,NML=quantum) CLOSE(In_Unit) DO iarran=1,narran nlegndre(iarran)=jmax(0,iarran) DO Vib=MinVib(iarran),MaxVib(iarran) jmax(Vib,iarran)=jmax(0,iarran) ENDDO ENDDO WRITE(Out_Unit,NML=quantum) OPEN(Unit=In_Unit,File=InputDIR(1:LEN(TRIM(InputDIR)))//InputFile,status='old') READ(In_Unit,NML=gauss) noscil=maxvib+1 CLOSE(In_Unit) WRITE(Out_Unit,NML=gauss) IF(NTheta_Plot<=0 .OR. NTheta_Plot>mxang)THEN WRITE(Out_Unit,*)'Error NTheta_Plot =', NTheta_Plot STOP 'error' ENDIF IF(NChi_Plot<=0 .OR. NChi_Plot>mxang)THEN WRITE(Out_Unit,*)'Error NChi_Plot =', NChi_Plot STOP 'error' ENDIF dtheta = halfpi/(NTheta_Plot-1) dchi= 4.d0*halfpi/(NChi_Plot-1) DO ichi = 1, NChi_Plot chi = dchi*(ichi-1) ChiOnly(ichi)=chi DO itheta = 1, NTheta_Plot Theta=dtheta*(itheta-1) ThetaOnly(Itheta)=Theta Theta_Plot(itheta+(ichi-1)*NTheta_Plot) = Theta Chi_Plot(itheta+(ichi-1)*NTheta_Plot) = chi ENDDO ENDDO CALL readbas(rho, NFull_States, MegaVal, nang, Phi_ABM, ThetaVal, ChiVals, weight, vect, NABM_Sfunc, eigen(1,MegaVal), NSYM_States, ithrho) IF(NFull_States>mxbasis)THEN WRITE(Out_Unit,*)'mxbasis<NFull_States:', mxbasis, NFull_States WRITE(Msg_Unit,*)'mxbasis<NFull_States:', mxbasis, NFull_States STOP "SFPlot_ABM" ENDIF WRITE(Out_Unit,*)'nang,mxang,NFull_States,mxbasis' WRITE(Out_Unit,*)nang, mxang, NFull_States, mxbasis CALL bascal2(nthetaChi_Plot, Theta_Plot, Chi_Plot, NFull_States, phinow, pn, dpndx, oscil, doscdz, parityfn, bfaci) ! -------------------------------------------------------------- ! DO similarity transform to get overlaps of surface functions. ! snew=s*vect contains overlaps of surface fcns. ! ----------------------------------------------------------------- IF(primitive)THEN WRITE(Out_Unit,*)'Primitive basis functions being plotted' ELSE WRITE(Out_Unit,*)'Constructing Actual Basis Functions' !CALL sgemm_junk(nthetaChi_Plot,NSym_States,NABM_Sfunc,phinow,nthetaChi_Plot,vect,NSym_States,sfun,NSym_States,0,1) CALL DGEMM('N','N',nthetaChi_Plot,NABM_Sfunc,NSym_States,1.d0,phinow,nthetaChi_Plot,vect,NSym_States,0.d0,sfun,nthetaChi_Plot) CALL vcopy(nthetaChi_Plot*NABM_Sfunc,phinow,sfun) ENDIF WRITE(Out_Unit,*)'ithmode=', ithmode IF(ithmode<=0)THEN WRITE(Out_Unit,*)'ithmode must be greater than 0' WRITE(Msg_Unit,*)'ithmode must be greater than 0' STOP "SFPlot_ABM" ELSEIF(primitive .AND. ithmode>NFull_States)THEN WRITE(Out_Unit,*)'ithmode must be less than NFull_States' WRITE(Msg_Unit,*)'ithmode must be less than NFull_States' STOP "SFPlot_ABM" ELSEIF(.NOT.primitive .AND. ithmode>NABM_Sfunc)THEN WRITE(Out_Unit,*)'ithmode must be less than NABM_Sfunc' WRITE(msg_Unit,*)'ithmode must be less than NABM_Sfunc' STOP "SFPlot_ABM" ENDIF RhoName =blank IF(ithrho<10)THEN WRITE(RhoName,'(i1)') ithrho ELSEIF(ithrho<100)THEN WRITE(RhoName,'(i2)') ithrho ELSEIF(ithrho<1000)THEN WRITE(RhoName,'(i3)') ithrho ELSEIF(ithrho<10000)THEN WRITE(RhoName,'(i4)') ithrho ELSEIF(ithrho<100000)THEN WRITE(RhoName,'(i5)') ithrho ELSE WRITE(Msg_Unit,*)'ithrho to large=',ithrho STOP 'SFPlot_ABM' ENDIF MegaName =blank IF(MegaVal<10)THEN WRITE(MegaName,'(i1)') MegaVal ELSEIF(MegaVal<100)THEN WRITE(MegaName,'(i2)') MegaVal ELSEIF(MegaVal<1000)THEN WRITE(MegaName,'(i3)') MegaVal ELSEIF(MegaVal<10000)THEN WRITE(MegaName,'(i4)') MegaVal ELSEIF(MegaVal<100000)THEN WRITE(MegaName,'(i5)') MegaVal ELSE WRITE(Msg_Unit,*)'MegaVal to large=',MegaVal STOP 'SFPlot_ABM' ENDIF OPEN(Unit=Plt_Unit,File=OutDIR(1:LEN(TRIM(OutDIR)))//'SfunctionPlots/Sfunc_ABM_'//Rhoname(1:LEN(TRIM(RhoName)))//'_'//Meganame(1:LEN(TRIM(MegaName)))//'.txt') valmax = 0.D0 SignAtMax = 0.D0 DO ichi = 1, NChi_Plot DO itheta = 1, NTheta_Plot phi = phinow(itheta+(ichi-1)*NTheta_Plot,ithmode) IF(ValMax<ABS(Phi))THEN ValMax=ABS(Phi) SignATMax=SIGN(1.d0,Phi) ENDIF ENDDO ENDDO DO ichi = 1, NChi_Plot DO itheta = 1, NTheta_Plot theta = thetaonly(itheta) chi = chionly(ichi) phi = phinow(itheta+(ichi-1)*NTheta_Plot,ithmode) CALL PMap(1,1,Zero,theta,chi,u,v) SP_PSI=phi*SignAtMax/ValMax IF(ABS(SP_PSI)<EPSILON(1.0_SP_Kind))SP_PSI=0.0 WRITE(Plt_Unit,'(3ES12.4)')u,v,SP_PSI Phi_Plot(itheta+(ichi-1)*NTheta_Plot) = phi*SignAtMax/ValMax ENDDO ENDDO CLOSE(Plt_Unit) WRITE(Out_Unit,99004)rho WRITE(Out_Unit,99004)1.1*valmax*cos(chif1(1)), 1.1*valmax*sin(chif1(1)), zero WRITE(Out_Unit,99004)1.1*valmax*cos(chif1(2)), 1.1*valmax*sin(chif1(2)), zero WRITE(Out_Unit,99004)1.1*valmax*cos(chif1(3)), 1.1*valmax*sin(chif1(3)), zero WRITE(Spin_Unit,99004)rho WRITE(Spin_Unit,99004)1.1*valmax*cos(chif1(1)), 1.1*valmax*sin(chif1(1)), zero WRITE(Spin_Unit,99004)1.1*valmax*cos(chif1(2)), 1.1*valmax*sin(chif1(2)), zero WRITE(Spin_Unit,99004)1.1*valmax*cos(chif1(3)), 1.1*valmax*sin(chif1(3)), zero DO ienergy = 1, 10 phi = energy(ienergy) IF(phi>phi_s)phi = phi_s*(1.0+c*log(phi/phi_s)) phi = (1.0+b*phi)*a WRITE(Spin_Unit,99004)phi, energy(ienergy) ENDDO CALL VecDef(thetaonly,NTheta_Plot,chionly,NChi_Plot,Phi_Plot,valmax) DEALLOCATE(thetaval) DEALLOCATE(chivals) DEALLOCATE(Phi_Plot) DEALLOCATE(Theta_Plot) DEALLOCATE(Chi_Plot) DEALLOCATE(weight) DEALLOCATE(phinow) DEALLOCATE(sfun) DEALLOCATE(parityfn) DEALLOCATE(bfaci) DEALLOCATE(pn) DEALLOCATE(dpndx) DEALLOCATE(oscil) DEALLOCATE(doscdz) DEALLOCATE(eigen) DEALLOCATE(vect) CLOSE(Out_Unit) OPEN(Unit=Out_Unit,File=OutDIR(1:LEN(TRIM(OutDIR)))//'Output/BasisForScatt.txt',Form='FORMATTED',ACCESS='APPEND') WRITE(Out_Unit,*)'Completed SFPlot_ABM_Main' WRITE(Out_Unit,*) RETURN 99001 FORMAT(a25) 99002 FORMAT(a25) 99003 FORMAT(1x,a70) 99004 FORMAT(1x,3E16.7) ENDSUBROUTINE SFPLOT_ABM_MAIN