print also energy difference relative to initial state

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.

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: ', &

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

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'