From 00d2d7eca30ee89cfb1a30c83f83cd80e0e753f9 Mon Sep 17 00:00:00 2001
From: Paul-Gerhard Reinhard <mpt218@homework.(none)>
Date: Mon, 9 May 2016 15:59:54 +0200
Subject: [PATCH] dipole excitations, some analyzing routines
---
Code/Changes.txt | 1 -
Code/Changes_private.txt | 10 ++
Code/external.f90 | 22 +++--
Code/util_PGR.f90 | 192 +++++++++++++++++++++++++++++++++++++++
4 files changed, 216 insertions(+), 9 deletions(-)
create mode 100644 Code/Changes_private.txt
create mode 100644 Code/util_PGR.f90
diff --git a/Code/Changes.txt b/Code/Changes.txt
index 7329861..9f68aba 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 0000000..65a310c
--- /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 1eac1c0..e2ff89f 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 0000000..2323e8d
--- /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
--
GitLab