From 497865d88788af9f9838c4c0a6cf2f7a2015fd8b Mon Sep 17 00:00:00 2001
From: Paul-Gerhard Reinhard <mpt218@homework.(none)>
Date: Fri, 20 May 2016 14:27:03 +0200
Subject: [PATCH] energy from external field, spectral analysis in 'Utils'

---
 Code/Makefile.ifort                 |  81 ++++++++++++++++++++++
 Code/Makefile.openmp                |   6 +-
 Code/dynamic.f90                    |  13 ++--
 Code/energies.f90                   |   1 +
 Code/external.f90                   |  16 ++++-
 Utils/Spectra/spectral_analysis.f90 | 103 ++++++++++++++++++++++++++++
 6 files changed, 212 insertions(+), 8 deletions(-)
 create mode 100644 Code/Makefile.ifort
 create mode 100644 Utils/Spectra/spectral_analysis.f90

diff --git a/Code/Makefile.ifort b/Code/Makefile.ifort
new file mode 100644
index 0000000..e3f4b01
--- /dev/null
+++ b/Code/Makefile.ifort
@@ -0,0 +1,81 @@
+COMPILER = ifort
+COMPILERFLAGS  = -axsse4.2 -msse4.2 -O3 -ip -no-prec-div -align all -mcmodel=large -shared-intel
+# -mkl=sequential
+LINKERFLAGS    =  $(COMPILERFLAGS)
+
+SOURCE_TRAILER = .f90
+EXEC           = sky3d.seq
+OBJS           = params.o grids.o levels.o fourier.o forces.o sequential.o \
+	inout.o coulomb.o trivial.o densities.o fragments.o twobody.o \
+	energies.o static.o meanfield.o dynamic.o pairs.o moment.o \
+	main3d.o user.o external.o abso_bc.o
+
+SHELL = /bin/sh
+
+#---------------------------------------------------------	
+all:$(EXEC) 
+$(EXEC):$(OBJS) 	
+	$(COMPILER) $(LINKERFLAGS) -o $@  $(OBJS)  -lfftw3 -llapack
+#
+
+user.o : user.f90 params.o grids.o levels.o
+
+external.o : external.f90 params.o sequential.o grids.o levels.o meanfield.o
+
+coulomb.o : coulomb.f90 params.o grids.o densities.o
+
+inout.o : inout.f90 params.o sequential.o grids.o moment.o densities.o \
+	levels.o  coulomb.o forces.o meanfield.o forces.o
+
+trivial.o : trivial.f90 params.o grids.o
+
+forces.o : forces.f90 forces.data params.o 
+
+fourier.o : fourier.f90 params.o grids.o
+
+params.o : params.f90
+
+moment.o : moment.f90 params.o grids.o densities.o
+
+sequential.o: sequential.f90 levels.o params.o grids.o densities.o
+
+levels.o : levels.f90 params.o grids.o forces.o fourier.o trivial.o
+
+grids.o : grids.f90 params.o forces.o
+
+densities.o : densities.f90 params.o grids.o levels.o trivial.o 
+
+fragments.o : fragments.f90 params.o grids.o forces.o levels.o twobody.o \
+	sequential.o 
+
+abso_bc.o : abso_bc.f90 params.o grids.o forces.o levels.o meanfield.o \
+	fragments.o sequential.o  Makefile
+
+twobody.o : twobody.f90 params.o grids.o densities.o moment.o forces.o
+
+energies.o :   energies.f90 params.o forces.o densities.o levels.o grids.o \
+	trivial.o pairs.o
+
+meanfield.o : meanfield.f90 params.o densities.o forces.o grids.o coulomb.o \
+	trivial.o levels.o
+
+pairs.o : pairs.f90 forces.o grids.o levels.o densities.o
+
+static.o : static.f90 params.o densities.o grids.o \
+	levels.o moment.o energies.o inout.o pairs.o meanfield.o
+
+dynamic.o : dynamic.f90 params.o densities.o grids.o trivial.o \
+	levels.o moment.o energies.o inout.o meanfield.o sequential.o \
+	twobody.o external.o abso_bc.o
+
+main3d.o : main3d.f90 params.o fourier.o forces.o densities.o meanfield.o levels.o grids.o fragments.o sequential.o dynamic.o static.o coulomb.o user.o
+
+.SUFFIXES:  .f90 .f .o 
+
+%.o : %.mod
+
+.f90.o:	
+	$(COMPILER) $(COMPILERFLAGS)  -c $<
+
+.f.o:
+	$(COMPILER) $(COMPILERFLAGS)  -c $<
diff --git a/Code/Makefile.openmp b/Code/Makefile.openmp
index e3e7209..8d2f8eb 100644
--- a/Code/Makefile.openmp
+++ b/Code/Makefile.openmp
@@ -1,13 +1,13 @@
 COMPILER = gfortran
-COMPILERFLAGS  =  -O3 -fopenmp
-LINKERFLAGS    =  -O3 -fopenmp
+COMPILERFLAGS  =  -O3 -fopenmp -msse4.2 -mfpmath=sse -ffast-math -finline-functions -funroll-loops
+LINKERFLAGS    =  $(COMPILERFLAGS)
 
 SOURCE_TRAILER = .f90
 EXEC           = sky3d.omp
 OBJS           = params.o grids.o levels.o fourier.o forces.o sequential.o \
 	inout.o coulomb.o trivial.o densities.o fragments.o twobody.o \
 	energies.o static.o meanfield.o dynamic.o pairs.o moment.o \
-	main3d.o user.o external.o
+	main3d.o user.o external.o abso_bc.o
 
 SHELL = /bin/sh
 
diff --git a/Code/dynamic.f90 b/Code/dynamic.f90
index 2897b93..b6d0aaa 100644
--- a/Code/dynamic.f90
+++ b/Code/dynamic.f90
@@ -66,7 +66,7 @@ CONTAINS
        ! Initialize *.res files
        CALL start_protocol(energiesfile, &
             '#    Time    N(n)    N(p)       E(sum)        E(integ)      Ekin &
-            &      Ecoll(n)     Ecoll(p)')
+            &      Ecoll(n)     Ecoll(p)  E_ext')
        CALL start_protocol(diffenergiesfile, &
             '#    Time N_n(0)-N_n(t) N_p(0)-N_p(t) diff-E(sum)  diff-E(integ)    diff-Ekin &
             &  diff-Ecoll(n)  diff-Ecoll(p)')
@@ -156,7 +156,7 @@ CONTAINS
                ps4,rho,tau,current,sdens,sodens)  
           psi(:,:,:,:,nst)=ps4
        ENDDO
-       !$OMP END PARALLEL DOx
+       !$OMP END PARALLEL DO
        ! sum up over nodes
        IF(tmpi) CALL collect_densities
        IF(nabsorb > 0) CALL absbc(nabsorb,iter,nt,time)
@@ -287,8 +287,13 @@ CONTAINS
        WRITE(scratch,'(1x,i5,9F10.4)') iter,orbital,spin,total_angmom 
        CLOSE(unit=scratch)
        OPEN(unit=scratch,file=energiesfile,POSITION='APPEND')  
-       WRITE(scratch,'(F10.2,2F8.3,2F15.7,3(1PG13.5))') &
-            time,pnr,ehf,ehfint,tke,ecoll
+       IF(ipulse==0) THEN  
+         WRITE(scratch,'(F10.2,2F8.3,2F15.7,3(1PG13.5))') &
+              time,pnr,ehf,ehfint,tke,ecoll
+       ELSE
+         WRITE(scratch,'(F10.2,2F8.3,2F15.7,4(1PG13.5))') &
+              time,pnr,ehf,ehfint,tke,ecoll,e_extern
+       END IF
        CLOSE(unit=scratch)
        OPEN(unit=scratch,file=diffenergiesfile,POSITION='APPEND')  
        WRITE(scratch,'(F10.2,9(1pg13.5))') &
diff --git a/Code/energies.f90 b/Code/energies.f90
index 404fadd..ba913b5 100644
--- a/Code/energies.f90
+++ b/Code/energies.f90
@@ -25,6 +25,7 @@ MODULE Energies
   REAL(db) :: ehf       ! Hartree-Fock energy from s.p. levels
   REAL(db) :: ehfprev
   REAL(db) :: e3corr    ! rearrangement energy
+  REAL(db) :: e_extern=0D0  ! energy from external field (accumulator)
   REAL(db) :: orbital(3),spin(3),total_angmom(3)
 CONTAINS
   !***************************************************************************
diff --git a/Code/external.f90 b/Code/external.f90
index e2ff89f..e245a63 100644
--- a/Code/external.f90
+++ b/Code/external.f90
@@ -4,8 +4,11 @@ MODULE External
   USE Grids, ONLY: nx,ny,nz,x,y,z,dx,dy,dz,wxyz
   USE Levels, ONLY: nstloc,isospin,charge_number,mass_number
   USE MEANFIELD, ONLY: upot
+  USE Densities, ONLY: rho
+  USE Energies, ONLY: e_extern
   IMPLICIT NONE  
-  INTEGER,PRIVATE :: isoext=0,ipulse=0
+  INTEGER :: ipulse=0
+  INTEGER,PRIVATE :: isoext=0
   REAL(db),PRIVATE,ALLOCATABLE,DIMENSION(:,:,:,:) :: extfield
   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
@@ -83,6 +86,10 @@ CONTAINS
     REAL(db) :: time
     INTENT(IN) :: time
     REAL(db) :: time_factor
+!   variables for computing energy absorbed from external field
+    REAL(db),SAVE :: time_factor_old=0D0
+    REAL(db),SAVE :: upot_ext_int,upot_ext_int_old=0D0
+!
     IF(ipulse==1) THEN  
        time_factor=EXP(-((time-tau0)/taut) **2)  
     ELSE  
@@ -98,6 +105,13 @@ CONTAINS
        time_factor=COS(omega*(time-tau0))*time_factor
     ENDIF
     upot=time_factor*extfield + upot
+
+    upot_ext_int = wxyz*SUM(extfield*rho)
+    e_extern=e_extern+(upot_ext_int+upot_ext_int_old)*&
+                      (time_factor-time_factor_old)
+    upot_ext_int_old = upot_ext_int
+    time_factor_old = time_factor
+
   END SUBROUTINE extfld
   !***********************************************************************
   SUBROUTINE extboost(noboostflag)
diff --git a/Utils/Spectra/spectral_analysis.f90 b/Utils/Spectra/spectral_analysis.f90
new file mode 100644
index 0000000..fe16f97
--- /dev/null
+++ b/Utils/Spectra/spectral_analysis.f90
@@ -0,0 +1,103 @@
+PROGRAM spectral_analysis
+  !
+  ! This program computes the spectral distribution of quadrupole 
+  ! oscillations produced as output on file 'quadrupoles.res'. 
+  ! Isoscalar and isovector quadrupole moments are composed from the
+  ! information on x**2, y**2, and z**2 in 'quadrupoles.res'.
+  ! A limiting time profile of type COS**n is multiplied to the signal
+  ! to guarantee that the signal approaches zero at the end of the
+  ! time interval (called time filtering, or windowing).
+  ! It is assumed that the quadrupole motion was initialized by an
+  ! instantaneous isoscalar quadrupole boost. The spectral strength is 
+  ! then simply the imaginary part of the Fourier transform of the 
+  ! isoscalar signal. The code also produces the Fourier transform
+  ! of the isovector signal. This has to be taken with care because
+  ! the excitation operator and the analzing operator are not the same
+  ! (what they should be for spectral analysis).
+  !
+  IMPLICIT NONE
+  INTEGER,PARAMETER :: db=SELECTED_REAL_KIND(12,100)
+  REAL(db),PARAMETER :: PI=3.14159265358979d0
+  REAL(db),PARAMETER :: hbc=197.32D0        ! to convert 1/fm to MeV
+  INTEGER,PARAMETER :: nfilter=4       ! order of damping COS**nfilter
+  !
+  REAL(db), ALLOCATABLE :: qext(:)
+  REAL(db), ALLOCATABLE :: qfilter(:)
+  REAL(db) :: time,timax,timin,deltime,timold,filter
+  REAL(db) :: delomega,omega
+  COMPLEX(db) :: cfac,cprod,cssum,cvsum,ctsum
+  INTEGER :: ndim
+  INTEGER :: nt,iomega
+
+
+  !
+  ! A first quick reading to check the length and integrity of 
+  ! the data file.
+  !
+  OPEN(1,file='extfield.res')
+  READ(1,'(1x)')
+  READ(1,'(1x)')
+  DO nt=1,999999
+     READ(1,*,ERR=19,END=19) timax
+     IF(nt==1) timin = timax
+     IF(nt==2) deltime=timax-timold
+     IF(nt>2 .AND. ABS(timax-timold-deltime)>1D-10) THEN
+        WRITE(*,*) 'problem at line ',nt+1
+        STOP "time steps not equidistant in 'extfield.res'"
+     END IF
+     ndim = nt
+     timold = timax
+  END DO
+19 CONTINUE
+  REWIND(1)
+  WRITE(*,'(a,i10,f12.2)') 'nr. of time steps, max.time=',ndim,timax
+  !  ALLOCATE(x2p(ndim),y2p(ndim),z2p(ndim),x2n(ndim),y2n(ndim),z2n(ndim))
+  ALLOCATE(qext(ndim),qfilter(ndim))
+
+  !
+  ! Read data and compose quadrupole moments (with time filter).
+  ! The coupling to isovector moment assumes an N=Z nucleus.
+  !
+  OPEN(2,file='extfield_filt.res')
+  WRITE(2,'(a)') '# external field signal in time domain:'
+  WRITE(2,'(a)') '# time[fm/c]     raw      filtered '
+  READ(1,'(1x)')
+  READ(1,'(1x)')
+  DO nt=1,ndim
+     READ(1,*,ERR=29,END=29) timax,qext(nt)
+     filter = COS((nt-1)*PI/(2D0*(ndim-1)))**nfilter    
+     qfilter(nt) = qext(nt)*filter
+     WRITE(2,'(f10.2,4(1pg13.5))') (nt-1)*deltime,qfilter(nt)
+  END DO
+  CLOSE(1)
+  CLOSE(2)
+  WRITE(*,*) "filtered signal written on 'extfield_filt.res''"
+
+  !
+  ! Perform Fourier transformation and print.
+  !
+  OPEN(2,file='extfield_spectrum.res')
+  WRITE(2,'(a)') '# spectral distributions of external-field strenght and power:'
+  WRITE(2,'(a)') '# omega[MeV]  REAL(FT signal)  strength     power'
+  deltime   = (timax-timin)/(ndim-1)
+  delomega  = (PI+PI)/(ndim*deltime)
+  DO iomega=1,ndim/4
+     omega  = iomega*delomega-delomega
+     cfac   = CEXP(CMPLX(0D0,omega*deltime))
+     cprod  = CMPLX(1D0,0D0)/cfac
+     cssum  = CMPLX(0D0,0D0)
+     DO nt=1,ndim
+        cprod = cprod*cfac
+        cssum = cprod*CMPLX(-qfilter(nt),0D0) + cssum
+     END DO
+     cssum = cssum*deltime
+     WRITE(2,'(f11.5,6(1pg13.5))') omega*hbc,cssum,ABS(cssum)**2
+  END DO
+  CLOSE(2)
+  WRITE(*,*) "spectra written on 'extfield_spectrum.res''"
+
+  STOP "regular end"
+
+29 STOP "error in reading data"
+
+END PROGRAM spectral_analysis
-- 
GitLab