dipole excitations, some analyzing routines

Code/Changes.txt

@@ -44,5 +44,4 @@ Step 2.0:
 
 
 
-
   

Code/Changes_private.txt

+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.
+  

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  

Code/util_PGR.f90

+! 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