diff --git a/Code/Changes.txt b/Code/Changes.txt index 732986126cbfea6b2ffc2ecbf315f8b79861fbdb..9f68abaac329cc91cf7819eec58bd3c45b694adf 100644 --- a/Code/Changes.txt +++ b/Code/Changes.txt @@ -44,5 +44,4 @@ Step 2.0: - diff --git a/Code/Changes_private.txt b/Code/Changes_private.txt new file mode 100644 index 0000000000000000000000000000000000000000..65a310cae4362040955fce0266e51e9565b7f8ef --- /dev/null +++ b/Code/Changes_private.txt @@ -0,0 +1,10 @@ +This file complements 'Changes.txt' with a protocol of changes +for the private branch 'master' only: + + +- Dipole excitations inserted in 'external.f90'. + +- Some analyzing routines (spin dipole, momentum space properties) + collected in 'user_PGR.f90'. These need explicit calls and coding + to be activated. + diff --git a/Code/external.f90 b/Code/external.f90 index 1eac1c0c3a2faa0f2e2711ae9161741a07d4130c..e2ff89f3ae915e1d29a5da4c2e40c31c4c3c83e9 100644 --- a/Code/external.f90 +++ b/Code/external.f90 @@ -7,15 +7,15 @@ MODULE External IMPLICIT NONE INTEGER,PRIVATE :: isoext=0,ipulse=0 REAL(db),PRIVATE,ALLOCATABLE,DIMENSION(:,:,:,:) :: extfield - REAL(db),PRIVATE :: amplq0=0.D0,radext=100.D0,widext=1.D0, & - tau0,taut,omega=0.D0 + REAL(db) :: amplq0=0.D0,amplx=0.D0,amply=0.D0,amplz=0.D0,& + radext=100.D0,widext=1.D0,tau0,taut,omega=0.D0 LOGICAL,PRIVATE :: textfield_periodic=.true. SAVE CONTAINS !*********************************************************************** SUBROUTINE getin_external - NAMELIST/extern/ amplq0,radext,widext,isoext,ipulse,omega,tau0,taut, & - textfield_periodic + NAMELIST/extern/ amplq0,amplx,amply,amplz,radext,widext,isoext,& + ipulse,omega,tau0,taut,textfield_periodic READ(5,extern) END SUBROUTINE getin_external SUBROUTINE init_external @@ -27,7 +27,9 @@ CONTAINS ' External field: called with invalid pulse type' IF(wflag) THEN WRITE(*,*) "***** Parameters of external field *****" - WRITE(*,"(a,e12.4)") " Amplitude of axial quad. =",amplq0 + WRITE(*,"(a,1pg12.4)") " Amplitude of axial quadripole =",amplq0 + WRITE(*,"(a,3(1pg12.4))") " Amplitudes Cartesian dipoles =",& + amplx,amply,amplz WRITE(*,"(2(A,F10.4),A)") " Radial damping: radius ",radext, & ' fm, width ',widext,' fm' WRITE(*,"(2(A,I2),2A)") " Isospin of excitation:",isoext, & @@ -61,10 +63,14 @@ CONTAINS DO ix=1,nx IF(textfield_periodic) THEN ! strictly periodic version facr=amplq0 *(2.D0*SIN(z(iz)*PI/zlim)**2 & - -SIN(x(ix)*PI/xlim)**2-SIN(y(iy)*PI/ylim)**2) + -SIN(x(ix)*PI/xlim)**2-SIN(y(iy)*PI/ylim)**2) & + +amplx*SIN(x(ix)*2D0*PI/xlim) & + +amply*SIN(y(iy)*2D0*PI/ylim) & + +amplz*SIN(z(iz)*2D0*PI/zlim) ELSE ! damped version - facr=amplq0 *(2.D0*z(iz)**2-x(ix)**2-y(iy)**2) & - /(1.0D0+EXP((SQRT(x(ix)**2+y(iy)**2+z(iz)**2)-radext)/widext)) + facr=(amplq0 *(2.D0*z(iz)**2-x(ix)**2-y(iy)**2) & + +amplx*x(ix)+amply*x(iy)+amplz*x(iz)) & + /(1.0D0+EXP((SQRT(x(ix)**2+y(iy)**2+z(iz)**2)-radext)/widext)) END IF extfield(ix,iy,iz,1)=facr*facn extfield(ix,iy,iz,2)=facr*facp diff --git a/Code/util_PGR.f90 b/Code/util_PGR.f90 new file mode 100644 index 0000000000000000000000000000000000000000..2323e8d0e3d0c9f66fe74f255af4c3feef3d9bf8 --- /dev/null +++ b/Code/util_PGR.f90 @@ -0,0 +1,192 @@ +! This is a collection of subroutines for analysis during +! TDHF computation. Each subroutein can be activated by +! moving it to 'user.f90' and inserting a CALL at the +! appropriate place. + + + +!*************************************************************************** +! This routine computes the expectation value of the spin-dipole +! and writes the result to the file 'spindipole.res'. +! The routine can be placed immediately in 'util.f90'. +! + SUBROUTINE spindipole() + USE Grids, ONLY: wxyz,x,y,z + USE Params, ONLY: time + + REAL(db) :: spindip(3,3) + INTEGER :: ix,iy,iz + LOGICAL,SAVE :: tinit=.TRUE. + + IF(tinit) THEN + OPEN(unit=scratch,file='spindipole.res') + WRITE(scratch,'(a)') & + '# time r_1*sigma_1 r_2*sigma_1 r_3*sigma_1 r_1*sigma_2 ....' + CLOSE(unit=scratch) + tinit=.FALSE. + END IF + + spindip = 0D0 + DO iz=1,nz + DO iy=1,ny + DO ix=1,nx + spindip(1,:) = spindip(1,:) + x(ix)*sdens(ix,iy,iz,:,1) + spindip(2,:) = spindip(2,:) + y(iy)*sdens(ix,iy,iz,:,1) + spindip(3,:) = spindip(3,:) + z(iz)*sdens(ix,iy,iz,:,1) + END DO + END DO + END DO + spindip = wxyz*spindip + + OPEN(unit=scratch,file='spindipole.res', POSITION='APPEND') + WRITE(scratch,'(f10.2,9(1pg13.5))') time,spindip + CLOSE(unit=scratch) + + END SUBROUTINE spindipole + + +! This routine tests the ortho-normality of the wavefunctions +! in the array 'psi'. +! The subroutine was part of a test version of 'levels.f90' +! and uses variables available in 'levels.f90'. Some new 'USE' +! statements may be necessary before usign it in 'user.f90'. +! + SUBROUTINE test_orthonorm(deviat_norm) + USE Trivial, ONLY: rpsnorm, overlap + REAL(db),INTENT(OUT) :: deviat_norm(2) + COMPLEX(db) :: oij + REAL(db) :: acc + INTEGER :: nst1,nst2,j,iq + DO iq=1,2 + acc = 0D0 + DO nst1=npmin(iq),npsi(iq) + DO nst2=nst1,npsi(iq) + oij=overlap(psi(:,:,:,:,nst1),psi(:,:,:,:,nst2)) + IF(nst1==nst2) THEN + acc = acc + (REAL(oij)-1D0)**2 + ELSE + acc = acc + 2D0*(REAL(oij)**2+AIMAG(oij)**2) + END IF + END DO + END DO + deviat_norm(iq) = SQRT(acc/(npsi(iq)-npmin(iq)+1)**2) + END DO + END SUBROUTINE test_orthonorm + + +! This routine computes the wavefunctions in momentum space +! and prints their absolute values at maxmimum momentum. This +! serves to check whether dynamics hits the upper bounds of +! momentum space. Results are written on unit 923 (unnamed). +! The subroutine was part of a test version of 'levels.f90' +! and uses variables available in 'levels.f90'. Some new 'USE' +! statements may be necessary before usign it in 'user.f90'. +! + SUBROUTINE pmomsmax(psin,timact) + COMPLEX(db), INTENT(IN) :: psin(:,:,:,:) + REAL(db), INTENT(IN) :: timact + COMPLEX(db), ALLOCATABLE:: ps1(:,:,:,:),ps2(:,:,:,:) + REAL(db) :: kfac + INTEGER :: iy + ALLOCATE(ps1(nx,ny,nz,2),ps2(nx,ny,nz,2)) + kfac=(PI+PI)/(dy*ny) + CALL dfftw_execute_dft(xforward,psin,ps1) + CALL dfftw_execute_dft(yforward,ps1,ps2) + CALL dfftw_execute_dft(zforward,ps2,ps1) + WRITE(921,'(f8.3,4(1pg13.5))') timact,SUM(ABS(ps1(nx/2,ny/2,nz/2,:))),& + SUM(ABS(ps1(nx/2,1,1,:))),SUM(ABS(ps1(1,ny/2,1,:))),& + SUM(ABS(ps1(1,1,nz/2,:))) + CALL flush(921) + DEALLOCATE(ps1,ps2) + END SUBROUTINE pmomsmax + + +! This routine computes the wavefunctions in momentum space and +! prints their absolute values along the bounds in momentum space +! in k_x, k_y, k_z direction. It serves to check on a broader +! basis whether and where dynamics hits the upper bounds of +! momentum space. Results are written on unit 922 (unnamed). +! The subroutine was part of a test version of 'levels.f90' +! and uses variables available in 'levels.f90'. Some new 'USE' +! statements may be necessary before usign it in 'user.f90'. +! + SUBROUTINE pmomscut(psin) + COMPLEX(db), INTENT(IN) :: psin(:,:,:,:) + COMPLEX(db), ALLOCATABLE:: ps1(:,:,:,:),ps2(:,:,:,:) + REAL(db) :: kfac + INTEGER :: iy,ix,iz + LOGICAL, SAVE :: tfirst=.true. + ! IF(.NOT.tfirst) RETURN + ALLOCATE(ps1(nx,ny,nz,2),ps2(nx,ny,nz,2)) + kfac=(PI+PI)/(dy*ny) + CALL dfftw_execute_dft(xforward,psin,ps1) + CALL dfftw_execute_dft(yforward,ps1,ps2) + CALL dfftw_execute_dft(zforward,ps2,ps1) + ! along x + WRITE(922,'(a)') '# along x' + DO ix=nx/2+1,nx + WRITE(922,'(f8.3,12(1pg13.5))') -(nx-ix+1)*kfac,& + (ABS(ps1(ix,ny/2,nz/2,:))),& + (ABS(ps1(ix,ny/2,1,:))),& + (ABS(ps1(ix,ny/2+1,1,:))),& + (ABS(ps1(ix,1,nz/2,:))),& + (ABS(ps1(ix,1,nz/2+1,:))),& + (ABS(ps1(ix,1,1,:))) + ENDDO + DO ix=1,nx/2 + WRITE(922,'(f8.3,12(1pg13.5))') (ix-1)*kfac,& + (ABS(ps1(ix,ny/2,nz/2,:))),& + (ABS(ps1(ix,ny/2,1,:))),& + (ABS(ps1(ix,ny/2+1,1,:))),& + (ABS(ps1(ix,1,nz/2,:))),& + (ABS(ps1(ix,1,nz/2+1,:))),& + (ABS(ps1(ix,1,1,:))) + ENDDO + WRITE(922,'(1x/1x)') + ! along y + WRITE(922,'(a)') '# along y' + DO iy=ny/2+1,ny + WRITE(922,'(f8.3,12(1pg13.5))') -(ny-iy+1)*kfac,& + (ABS(ps1(nx/2,iy,nz/2,:))),& + (ABS(ps1(nx/2,iy,1,:))),& + (ABS(ps1(nx/2,iy,nz/2+1,:))),& + (ABS(ps1(1,iy,nz/2,:))),& + (ABS(ps1(nx/2+1,iy,nz/2,:))),& + (ABS(ps1(1,iy,1,:))) + ENDDO + DO iy=1,ny/2 + WRITE(922,'(f8.3,12(1pg13.5))') (iy-1)*kfac,& + (ABS(ps1(nx/2,iy,nz/2,:))),& + (ABS(ps1(nx/2,iy,1,:))),& + (ABS(ps1(nx/2,iy,nz/2+1,:))),& + (ABS(ps1(1,iy,nz/2,:))),& + (ABS(ps1(nx/2+1,iy,nz/2,:))),& + (ABS(ps1(1,iy,1,:))) + ENDDO + WRITE(922,'(1x/1x)') + ! along z + WRITE(922,'(a)') '# along z' + DO iz=nz/2+1,nz + WRITE(922,'(f8.3,12(1pg13.5))') -(nz-iz+1)*kfac,& + (ABS(ps1(nx/2,ny/2,iz,:))),& + (ABS(ps1(1,ny/2,iz,:))),& + (ABS(ps1(nx/2+1,ny/2,iz,:))),& + (ABS(ps1(nx/2,1,iz,:))),& + (ABS(ps1(nx/2,ny/2+1,iz,:))),& + (ABS(ps1(1,1,iz,:))) + ENDDO + DO iz=1,nz/2 + WRITE(922,'(f8.3,12(1pg13.5))') (iz-1)*kfac,& + (ABS(ps1(nx/2,ny/2,iz,:))),& + (ABS(ps1(1,ny/2,iz,:))),& + (ABS(ps1(nx/2+1,ny/2,iz,:))),& + (ABS(ps1(nx/2,1,iz,:))),& + (ABS(ps1(nx/2,ny/2+1,iz,:))),& + (ABS(ps1(1,1,iz,:))) + ENDDO + ! + WRITE(922,'(1x/1x)') + CALL flush(922) + DEALLOCATE(ps1,ps2) + ! tfirst=.false. + END SUBROUTINE pmomscut