From f95b1046f3f95776af8ebd6b7f5aec5544009df4 Mon Sep 17 00:00:00 2001
From: Paul-Gerhard Reinhard <mpt218@homework.(none)>
Date: Wed, 11 May 2016 20:19:13 +0200
Subject: [PATCH] print also energy difference relative to initial state
---
Code/Changes_private.txt | 4 ++++
Code/dynamic.f90 | 21 +++++++++++++++++++--
Code/main3d.f90 | 3 ++-
Code/params.f90 | 1 +
4 files changed, 26 insertions(+), 3 deletions(-)
diff --git a/Code/Changes_private.txt b/Code/Changes_private.txt
index 65a310c..b52840a 100644
--- a/Code/Changes_private.txt
+++ b/Code/Changes_private.txt
@@ -8,3 +8,7 @@ for the private branch 'master' only:
collected in 'user_PGR.f90'. These need explicit calls and coding
to be activated.
+- absorbing boundary conditions implemented
+ (special obseravbles as PES and PAD yet to be done).
+
+- protocol of energy and particle-number difference to initial state.
diff --git a/Code/dynamic.f90 b/Code/dynamic.f90
index 915a7f2..2897b93 100644
--- a/Code/dynamic.f90
+++ b/Code/dynamic.f90
@@ -24,6 +24,8 @@ MODULE DYNAMIC
LOGICAL :: texternal=.FALSE. ! must be true if an external field is present
LOGICAL :: text_timedep ! true for time-dependent external field
REAL(db),PARAMETER :: esf=0.0D0
+ REAL(db) :: pnrold(2)=(/0D0,0D0/),ehfold=0D0,eintold=0D0,ekinold=0D0,&
+ ecollold(2)=0D0
CONTAINS
!*************************************************************************
SUBROUTINE getin_dynamic
@@ -49,7 +51,7 @@ CONTAINS
! Step 1: Preparation phase
IF(.NOT.trestart) THEN
iter=0
- time=0.0D0
+ time=0D0
! save wave functions
CALL write_wavefunctions
IF(wflag) WRITE(*,*) 'Wrote wave function file after initialization'
@@ -65,6 +67,9 @@ CONTAINS
CALL start_protocol(energiesfile, &
'# Time N(n) N(p) E(sum) E(integ) Ekin &
& Ecoll(n) Ecoll(p)')
+ 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)')
CALL start_protocol(monopolesfile, &
'# Time rms_n rms_p rms_tot rms_n-rms_p')
CALL start_protocol(dipolesfile, &
@@ -151,7 +156,7 @@ CONTAINS
ps4,rho,tau,current,sdens,sodens)
psi(:,:,:,:,nst)=ps4
ENDDO
- !$OMP END PARALLEL DO
+ !$OMP END PARALLEL DOx
! sum up over nodes
IF(tmpi) CALL collect_densities
IF(nabsorb > 0) CALL absbc(nabsorb,iter,nt,time)
@@ -271,6 +276,13 @@ CONTAINS
(current(:,:,:,1,iq)**2+current(:,:,:,2,iq)**2&
+current(:,:,:,3,iq)**2)/rho(:,:,:,iq) )
ENDDO
+ IF(time==0D0) THEN
+ pnrold=pnr
+ ehfold=ehf
+ eintold=ehfint
+ ekinold=tke
+ ecollold=ecoll
+ END IF
OPEN(unit=scratch,file=spinfile, POSITION='APPEND')
WRITE(scratch,'(1x,i5,9F10.4)') iter,orbital,spin,total_angmom
CLOSE(unit=scratch)
@@ -278,6 +290,11 @@ CONTAINS
WRITE(scratch,'(F10.2,2F8.3,2F15.7,3(1PG13.5))') &
time,pnr,ehf,ehfint,tke,ecoll
CLOSE(unit=scratch)
+ OPEN(unit=scratch,file=diffenergiesfile,POSITION='APPEND')
+ WRITE(scratch,'(F10.2,9(1pg13.5))') &
+ time,pnr-pnrold,ehf-ehfold,ehfint-eintold,&
+ tke-ekinold,ecoll-ecollold
+ CLOSE(unit=scratch)
WRITE(*,'(/A,I7,A,F12.4,A/2(A,F12.4),A)') &
' ***** Time step #',iter,' at time ',time,' fm/c *****', &
' Total energy: ',ehf,' MeV Total kinetic energy: ', &
diff --git a/Code/main3d.f90 b/Code/main3d.f90
index 15e0e9d..5a56f76 100644
--- a/Code/main3d.f90
+++ b/Code/main3d.f90
@@ -16,7 +16,8 @@ PROGRAM tdhf3d
INTEGER :: imode,nofsave
!***********************************************************************
NAMELIST /files/ wffile,converfile,monopolesfile,dipolesfile, &
- momentafile,energiesfile,quadrupolesfile,spinfile,extfieldfile
+ momentafile,energiesfile,quadrupolesfile,spinfile,extfieldfile,&
+ diffenergiesfile
NAMELIST /main/ tcoul,mprint,mplot,trestart, &
writeselect,write_isospin,mrest,imode,tfft,nof,r0
!********************************************************************
diff --git a/Code/params.f90 b/Code/params.f90
index 7a993e0..6a97da4 100644
--- a/Code/params.f90
+++ b/Code/params.f90
@@ -19,6 +19,7 @@ MODULE Params
CHARACTER(LEN=80) :: dipolesfile='dipoles.res'
CHARACTER(LEN=80) :: momentafile='momenta.res'
CHARACTER(LEN=80) :: energiesfile='energies.res'
+ CHARACTER(LEN=80) :: diffenergiesfile='diffenergies.res'
CHARACTER(LEN=80) :: quadrupolesfile='quadrupoles.res'
CHARACTER(LEN=80) :: spinfile='spin.res'
CHARACTER(LEN=80) :: extfieldfile='extfield.res'
--
GitLab