From 7d9659f00ff49e6ae374ba9594f32a32d7f8c816 Mon Sep 17 00:00:00 2001 From: Aaron Philip Date: Wed, 11 Mar 2026 13:18:32 -0400 Subject: [PATCH 1/8] Implemented absorbing conditions compatible with main branch of Sky3D --- Code/Makefile | 8 +++++--- Code/dynamic.f90 | 6 +++++- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/Code/Makefile b/Code/Makefile index 3082aed..56a8727 100644 --- a/Code/Makefile +++ b/Code/Makefile @@ -11,7 +11,7 @@ ifort ifort_seq ifort_omp : export COMPILER_SKY = ifort # # ----- define libraries----------------------------------------------------- # -seq omp debug seq_debug omp_debug ifort ifort_seq ifort_omp : export LIBS_SKY = -lfftw3 -llapack -lopenblas +seq omp debug seq_debug omp_debug ifort ifort_seq ifort_omp : export LIBS_SKY = -lfftw3 -lscalapack -lopenblas mpi mpi_debug mpi-omp : export LIBS_SKY = -lfftw3 -lscalapack -llapack -lopenblas apple apple_omp : export LIBS_SKY = -lfftw3 -framework Accelerate -L/opt/homebrew/lib # @@ -60,7 +60,7 @@ clean-exec: OBJS = params.o grids.o levels.o fourier.o forces.o $(PARALLEL_SKY) \ 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 ylm.o + main3d.o user.o external.o ylm.o abso_bc.o # LINKERFLAGS_SKY = $(COMPILERFLAGS_SKY) # @@ -124,7 +124,9 @@ static.o : static.f90 params.o densities.o grids.o \ dynamic.o : dynamic.f90 params.o densities.o grids.o trivial.o \ levels.o moment.o energies.o inout.o meanfield.o $(PARALLEL_SKY) \ - twobody.o external.o Makefile + twobody.o external.o abso_bc.o Makefile + +abso_bc.o : abso_bc.f90 params.o grids.o levels.o densities.o Makefile main3d.o : main3d.f90 params.o fourier.o forces.o densities.o meanfield.o \ levels.o grids.o fragments.o $(PARALLEL_SKY) dynamic.o static.o coulomb.o user.o \ diff --git a/Code/dynamic.f90 b/Code/dynamic.f90 index 80c3d80..e40306b 100644 --- a/Code/dynamic.f90 +++ b/Code/dynamic.f90 @@ -22,8 +22,10 @@ MODULE DYNAMIC USE Inout, ONLY: write_wavefunctions,write_densities, plot_density, & sp_properties,start_protocol USE External + USE abso_bc IMPLICIT NONE SAVE + INTEGER :: nabsorb=0 INTEGER :: nt !< the number of the final time step to be !! calculated. In case of a restart this is smaller than the total number of time steps. REAL(db) :: dt !< the physical time increment in units of fm/c. @@ -49,7 +51,7 @@ MODULE DYNAMIC !!parameters. !--------------------------------------------------------------------------- SUBROUTINE getin_dynamic - NAMELIST /dynamic/ nt,dt,mxpact,mrescm,rsep,texternal + NAMELIST /dynamic/ nt,dt,mxpact,mrescm,rsep,texternal,nabsorb READ(5,dynamic) IF(wflag) THEN WRITE(*,*) '***** Parameters for the dynamic calculation *****' @@ -58,6 +60,7 @@ SUBROUTINE getin_dynamic WRITE(*,'(A,F7.2,A)') ' The calculation stops at ',rsep, & ' fm fragment separation' WRITE(*,'(A,I3)') ' Power limit in operator expansion:',mxpact + WRITE(*,'(A,I3)') ' Number of absorbing points:',nabsorb ENDIF IF(texternal) CALL getin_external END SUBROUTINE getin_dynamic @@ -221,6 +224,7 @@ SUBROUTINE dynamichf !$OMP END PARALLEL DO IF(tmpi) CALL collect_densities ! calculate mean fields and external fields + IF(nabsorb > 0) CALL absbc(nabsorb,iter,nt,time) CALL skyrme IF(text_timedep) CALL extfld(0.D0) CALL tinfo From f8ca9ba6ddac80528460b34e67b17e732084226a Mon Sep 17 00:00:00 2001 From: Aaron Philip Date: Wed, 11 Mar 2026 13:30:18 -0400 Subject: [PATCH 2/8] Added missing abso_bc.f90 copied over from astro-parallel branch --- Code/abso_bc.f90 | 535 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 535 insertions(+) create mode 100644 Code/abso_bc.f90 diff --git a/Code/abso_bc.f90 b/Code/abso_bc.f90 new file mode 100644 index 0000000..f52baed --- /dev/null +++ b/Code/abso_bc.f90 @@ -0,0 +1,535 @@ +!------------------------------------------------------------------------------ +! MODULE: abso_bc +!------------------------------------------------------------------------------ +! DESCRIPTION: +!> @brief +!!Package containing absorbing boundary conditions and +!!subsequent analysis of observables from electron emission. +!!One may also consider to perform the cumulation of total +!!absorbed density within the absorbing loop. +!------------------------------------------------------------------------------ +MODULE abso_bc + USE Params + USE Grids, ONLY: dx,dy,dz,x,y,z,wxyz + USE Levels, ONLY: nstmax,psi,npsi,npmin,wocc,isospin,nprot,nneut + USE Meanfield + USE Parallel + IMPLICIT NONE + SAVE + !*********************************************************************** + ! * + ! General variables for the absorbing bounds * + ! * + !*********************************************************************** + INTEGER,PARAMETER :: nocc=1 !< degeneracy of states (outdated) + INTEGER :: ispherabso=0 !< switch to spherical absorbing bounds + INTEGER :: iangabso=0 + INTEGER :: ipes=0 + INTEGER :: ifabsorbital=0 !< switch to orbit-wise accumulation + INTEGER :: nangtheta=0 !< number of reference points for outgoing spectra + INTEGER :: nangphi=0 + INTEGER :: jescmaskorb=0 + INTEGER :: jescmask=10 + REAL(db) :: powabs=0.0175D0 !< power of absorbing mask + + REAL(db),PRIVATE :: xango + REAL(db),PRIVATE :: yango + REAL(db),PRIVATE :: zango + INTEGER,PRIVATE :: ix + INTEGER,PRIVATE :: iy + INTEGER,PRIVATE :: iz + INTEGER,PRIVATE :: nst + INTEGER,PRIVATE :: iq + INTEGER,PRIVATE :: nta + INTEGER,PRIVATE :: ita + INTEGER,PRIVATE :: nabso=0 ! @brief +!!Apply absorbing bounds, optionally accumulate absorbed density. +!--------------------------------------------------------------------------- + SUBROUTINE absbc(nabsorb,it,nt,time) + INTEGER,INTENT(in) :: it,nt,nabsorb + REAL(db),INTENT(in) :: time + REAL(db) :: weight + INTEGER :: iqa + + LOGICAL,SAVE :: firstcall = .TRUE. + REAL(db),ALLOCATABLE,DIMENSION(:,:,:,:) :: rhoabso_all + + !------------------------------------------------------------ + + nabso = nabsorb + IF(nabsorb <= 0) RETURN + + ! set calling parameters to internal common variables + ita = it + nta = nt + timea = time + + ! initializations + + IF (firstcall) THEN + IF(nabsorb>=NX/2.OR.nabsorb>=NY/2 .OR. nabsorb>=NZ/2) & + STOP 'ABSBC: problem with nabsorb' + CALL init_absbc() + firstcall = .FALSE. + END IF + + ! applying absorbing masks and accumulating absorption + + + + + ! apply mask function (and accumulate absorption per state) + ! and compute new density + + DO nst = 1,nstloc + iqa = isospin(globalindex(nst)) + weight = nocc*wocc(nst) + DO iz=1,NZ; DO iy=1,NY; DO ix=1,NX + IF(tgridabso(ix,iy,iz)) THEN + rhoabso(ix,iy,iz,iqa) = rhoabso(ix,iy,iz,iqa) + & + weight*spherloss(ix,iy,iz)*( & + abs2(psi(ix,iy,iz,1,nst))+abs2(psi(ix,iy,iz,2,nst)) ) + psi(ix,iy,iz,1,nst)=absomask(ix,iy,iz)*psi(ix,iy,iz,1,nst) + psi(ix,iy,iz,2,nst)=absomask(ix,iy,iz)*psi(ix,iy,iz,2,nst) + END IF + END DO;END DO;END DO + END DO + + IF(ifabsorbital==1) THEN + DO nst = 1,nstloc + DO iz=1,NZ; DO iy=1,NY; DO ix=1,NX + IF(tgridabso(ix,iy,iz)) THEN + rhoescmaskorb(ix,iy,iz,nst) = rhoescmaskorb(ix,iy,iz,nst) & + + spherloss(ix,iy,iz)*( & + abs2(psi(ix,iy,iz,1,nst))+abs2(psi(ix,iy,iz,2,nst)) ) + END IF + END DO;END DO;END DO + END DO + END IF + + + + ! print escaped nucleons + + IF(tmpi .AND. (MOD(ita,mprint)==1 .OR. MOD(ita,mrest)==0 .OR. & + MOD(ita,jescmask) == 0 .OR. ita == nta) ) THEN + ALLOCATE(rhoabso_all(NX,NY,NZ,2)) + CALL mpi_barrier (mpi_comm_world, mpi_ierror) + CALL mpi_allreduce(rhoabso,rhoabso_all,2*NX*NY*NZ, & + mpi_double_precision,mpi_sum, & + mpi_comm_world,mpi_ierror) + IF(MOD(ita,mprint)==1) CALL nescape(rhoabso_all) + CALL escmask(rhoabso_all) + ELSE + IF(MOD(ita,mprint)==1) CALL nescape(rhoabso) + CALL escmask(rhoabso) + END IF + + ! prints wavefunction on measuring points on file + + ! IF (ipes /= 0) CALL evalMP() ! not yet implemented + + IF(MOD(ita,mrest)==0) THEN + IF(tmpi) THEN + CALL wrtabso(ita,rhoabso_all) + ELSE + CALL wrtabso(ita,rhoabso) + END IF + WRITE(6,*) ' Augment restart file at it=',ita + ENDIF + IF(ALLOCATED(rhoabso_all)) DEALLOCATE(rhoabso_all) + + RETURN + END SUBROUTINE absbc +!--------------------------------------------------------------------------- +! DESCRIPTION: init_absbc +!> @brief +!!Initializes absorbing boundary conditions. +!--------------------------------------------------------------------------- + + SUBROUTINE init_absbc() + + NAMELIST /absobc/ ispherabso,iangabso,ipes,ifabsorbital, & + nangtheta,nangphi,jescmaskorb,jescmask,powabs + + !READ(5,absobc) + + ALLOCATE(tgridabso(NX,NY,NZ)) + ALLOCATE(absomask(NX,NY,NZ)) + ALLOCATE(spherloss(NX,NY,NZ)) + IF(ispherabso /= 0) THEN + CALL init_spherabso() + ELSE + CALL init_abso() + END IF + IF(ifabsorbital==1) THEN + ALLOCATE(rhoescmaskorb(NX,NY,NZ,nstloc)) + IF(.NOT.trestart) rhoescmaskorb = 0D0 + END IF + ALLOCATE(rhoabso(NX,NY,NZ,2)) + IF(.NOT.trestart) rhoabso = 0D0 + ! IF (ipes /= 0) CALL init_MP() ! not yet implemented + IF(trestart) THEN + CALL readabso() + WRITE(6,*) ' Reading abso-info from restart file' + ENDIF + + END SUBROUTINE init_absbc + + + + + REAL(db) FUNCTION dist_min() + + ! determines minimum radius for last non-absorbing point + + REAL(db) :: rmin + + xango=0D0 + yango=0D0 + zango=0D0 + + rmin = abs(x(NX-nabso)-xango) + rmin = min(rmin,abs(x(nabso+1)-xango)) + rmin = min(rmin,abs(y(NY-nabso)-yango)) + rmin = min(rmin,abs(y(nabso+1)-yango)) + rmin = min(rmin,abs(z(NZ-nabso)-zango)) + rmin = min(rmin,abs(z(nabso+1)-zango)) + + dist_min = rmin + + + RETURN + END FUNCTION dist_min +!--------------------------------------------------------------------------- +! DESCRIPTION: init_sphereabso +!> @brief +!!Initializes sphericalabsorbing boundary conditions. +!--------------------------------------------------------------------------- + SUBROUTINE init_spherabso() + + ! Initializes mask function for spherical boundary conditions + + REAL(db) :: dmin1,dmin2,bcrad,dmin12,dmin22,dist,dist2,cosact + + !------------------------------------------------------------ + + bcrad = nabso*dx + dmin1 = dist_min() + IF (DMIN1 < 0.) STOP 'Error in abso: dmin1<0' + dmin2 = dmin1+bcrad + dmin22=dmin2**2 + dmin12=DMIN1**2 + + DO iz=1,NZ; DO iy=1,NY; DO ix=1,NX + dist2 = x(ix)**2+y(iy)**2+z(iz)**2 + IF (dist2 <= dmin12) THEN + absomask(ix,iy,iz)=1D0 + tgridabso(ix,iy,iz)=.false. + ELSE IF (dist2 > dmin22) THEN + absomask(ix,iy,iz)=0D0 + tgridabso(ix,iy,iz)=.true. + ELSE + dist = MAX(1D-20,SQRT(dist2)) + cosact = COS((dist-dmin1)*0.5D0*pi/bcrad) + IF(cosact > 0D0) THEN + absomask(ix,iy,iz)= cosact**powabs + ELSE + absomask(ix,iy,iz)= 0D0 + END IF + tgridabso(ix,iy,iz)=.true. + END IF + spherloss(ix,iy,iz) = 1D0-absomask(ix,iy,iz)**2 + END DO;END DO;END DO + + RETURN + END SUBROUTINE init_spherabso +!--------------------------------------------------------------------------- +! DESCRIPTION: init_abso +!> @brief +!!Initializes normal absorbing boundary conditions. +!--------------------------------------------------------------------------- + SUBROUTINE init_abso() + USE Parallel + + ! Initializes mask for rectangular absorbing boundaries conditions + + REAL(db),ALLOCATABLE,DIMENSION(:) :: xmask,ymask,zmask + + LOGICAL :: wflagabs=.true. + + ALLOCATE(xmask(NX)) + ALLOCATE(ymask(NY)) + ALLOCATE(zmask(NZ)) + + ! prepare mask functions in each direction separately + + zmask = 1D0 + DO iz = 1,nabso + zmask(iz) = COS(pi*0.5D0*(nabso+1.0-iz)/nabso)**powabs + zmask(NZ+1-iz) = zmask(iz) + END DO + + ymask = 1D0 + DO iy = 1,nabso + ymask(iy) = COS(pi*0.5D0*(nabso+1.0-iy)/nabso)**powabs + ymask(NY+1-iy) = ymask(iy) + END DO + + xmask = 1D0 + DO ix = 1,nabso + xmask(ix) = COS(pi*0.5D0*(nabso+1.0D0-ix)/nabso)**powabs + xmask(NX+1-ix) = xmask(ix) + END DO + + IF(mpi_myproc == 0 .AND. wflagabs) THEN + WRITE(6,'(a)') ' ZMASK:' + WRITE(6,'(1x,5(1pg12.4))') zmask + WRITE(6,'(a)') ' YMASK:' + WRITE(6,'(1x,5(1pg12.4))') ymask + WRITE(6,'(a)') ' XMASK:' + WRITE(6,'(1x,5(1pg12.4))') xmask + END IF + + + ! compose to one mask function on all grid points + + FORALL(ix=1:NX,iy=1:NY,iz=1:NZ) + absomask(ix,iy,iz)=xmask(ix)*ymask(iy)*zmask(iz) + tgridabso(ix,iy,iz)= absomask(ix,iy,iz) < 0.999999999999D0 + spherloss(ix,iy,iz) = 1D0-absomask(ix,iy,iz)**2 + END FORALL + + RETURN + END SUBROUTINE init_abso + + + + !-----escmask--------------------------------------------------------- + + SUBROUTINE escmask(rhoabsoprint) + USE Parallel + + ! print collected information on escaping nucleons + + REAL(db),INTENT(IN),DIMENSION(:,:,:,:) :: rhoabsoprint + + CHARACTER(LEN=5) :: scratchline + CHARACTER(LEN=5) :: readstring + CHARACTER(LEN=4) :: str + + !-------------------------------------------------------------------- + + + + + IF(ifabsorbital==1) THEN + IF (MOD(ita,jescmaskorb) == 0 .OR. ita==nta) THEN + DO nst=1,nstmax + IF(nstmax > 9999) STOP 'ERROR: Too many states for escmaskOrb' + IF(mpi_myproc == node(nst)) THEN + WRITE(scratchline,'(i5)') 10000+nst + READ(scratchline,'(a)') readstring + str = readstring(2:5) + OPEN(588,STATUS='unknown', FILE='escmaskOrb.'//str) + WRITE(588,'(a,f12.4,a,i5,/a,3i5/a)') & + '# distribution of escaped nucleon at time=',timea,& + ' for state=',nst, & + '#nx,ny,nz:',NX,NY,NZ, & + '# x y z loss of density' + + DO iz=1,NZ; DO iy=1,NY; DO ix=1,NX + WRITE(588,'(3f12.4,1e17.7)') & + x(ix),y(iy),z(iz),rhoescmaskorb(ix,iy,iz,node(nst)) + END DO; END DO; END DO + + CLOSE(588) + END IF + END DO + END IF + END IF + + + IF(mpi_myproc == 0) THEN + IF (MOD(ita,jescmask) == 0 .OR. ita == nta) THEN + DO iq=1,2 + IF(iq == 1)THEN + OPEN(589,STATUS='unknown',FILE='escmask.1') + WRITE(589,'(a,f12.4/a,3i5/a)') & + '# total distribution of escaped protons at time=',timea, & + '#nx,ny,nz:',NX,NY,NZ, & + '# x y z loss of density' + ELSE + OPEN(589,STATUS='unknown',FILE='escmask.2') + WRITE(589,'(a,f12.4/a,3i5/a)') & + '# total distribution of escaped neutrons at time=',timea, & + '#nx,ny,nz:',NX,NY,NZ, & + '# x y z loss of density' + END IF + + DO iz=1,NZ; DO iy=1,NY; DO ix=1,NX + WRITE(589,'(3f12.4,1e17.7)') & + x(ix),y(iy),z(iz),rhoabsoprint(ix,iy,iz,iq) + END DO; END DO; END DO + + CLOSE(589) + END DO + END IF + END IF + + + RETURN + END SUBROUTINE escmask + + + + + SUBROUTINE nescape(rhoabsoprint) + + ! Compute and print total number of escaped nucleons + + REAL(db),INTENT(IN),DIMENSION(:,:,:,:) :: rhoabsoprint + + LOGICAL,SAVE :: firstcall = .TRUE. + LOGICAL :: texist + REAL(db) :: aprot,aneut + + !------------------------------------------------------------ + + IF(firstcall) THEN + INQUIRE(file='nescape.res',exist=texist) + IF(.NOT.texist) THEN + OPEN(unit=scratch,file='nescape.res') + WRITE(scratch,'(a)') & + '# number of escaped nucleons ', & + '# time lost protons lost neutrons ' + CLOSE(unit=scratch) + END IF + firstcall = .FALSE. + END IF + + aprot = 1D0*nprot !SUM(fcharg(1:2)) + aneut = 1D0*nneut !SUM(fmass(1:2))-aprot + IF(mpi_myproc == 0) THEN + OPEN(unit=scratch,file='nescape.res',POSITION='APPEND') + WRITE(scratch,'(1x,f10.4,4(1pg20.10))') timea, & + aprot-wxyz*sum(rho(:,:,:,1)),aneut-wxyz*sum(rho(:,:,:,2)), & + wxyz*sum(rhoabsoprint(:,:,:,1)),wxyz*sum(rhoabsoprint(:,:,:,2)) + CLOSE(unit=scratch) + END IF + + RETURN + END SUBROUTINE nescape + + !------------------------------------------------------------ + + + SUBROUTINE wrtabso(iteration,rhoabsoprint) + ! + USE Parallel, ONLY: mpi_myproc,mpi_nprocs,nstloc + IMPLICIT NONE + INTEGER,INTENT(IN) :: iteration + INTEGER :: nst,lenstr + CHARACTER(120) :: rsfp + REAL(db),INTENT(IN),DIMENSION(:,:,:,:) :: rhoabsoprint + ! + !*********************************************************************** + ! * + ! complement wffile by info on accumulated absorption * + ! * + !*********************************************************************** + ! + ! first write general information only with processor 0 + ! + IF(wffile=='none'.OR.wffile=='NONE') THEN + WRITE(6,*) " wffile='NONE' --> no file written " + RETURN + ENDIF + lenstr = LEN_TRIM(wffile) + IF(mpi_myproc==0) THEN + OPEN(UNIT=scratch2,FILE=wffile(1:lenstr)//'.a',STATUS='REPLACE', & + FORM='UNFORMATTED') + WRITE(scratch2) rhoabsoprint + ENDIF + ! now if we are running with one processor, also write the wave functions + ! onto that file; otherwise open new individualized file + ! + IF(ifabsorbital==1) THEN + IF(mpi_nprocs>1) THEN + IF(mpi_myproc/=0) CLOSE(UNIT=scratch2,STATUS='KEEP') + WRITE(rsfp,'(I3.3,''.'',A)') mpi_myproc,wffile(1:lenstr)//'.a' + OPEN(UNIT=scratch2,FILE=rsfp,STATUS='REPLACE', & + FORM='UNFORMATTED') + ENDIF + ! + DO nst = 1,nstloc + WRITE(scratch2) rhoescmaskOrb(:,:,:,nst) + ENDDO + END IF + CLOSE(UNIT=scratch2,STATUS='KEEP') + WRITE(6,*) 'Restart file augmented at iteration ',iteration + END SUBROUTINE wrtabso + + SUBROUTINE readabso() + ! + USE Parallel, ONLY: mpi_myproc,mpi_nprocs,nstloc + IMPLICIT NONE + INTEGER :: nst,lenstr + CHARACTER(120) :: rsfp + ! + !*********************************************************************** + ! * + ! complement wffile by info on accumulated absorption * + ! * + !*********************************************************************** + ! + ! first write general information only with processor 0 + ! + IF(wffile=='none'.OR.wffile=='NONE') THEN + WRITE(6,*) " wffile='NONE' --> no file read " + RETURN + ENDIF + lenstr = LEN_TRIM(wffile) + IF(mpi_myproc==0) THEN + OPEN(UNIT=scratch2,FILE=wffile(1:lenstr)//'.a',STATUS='OLD', & + FORM='UNFORMATTED') + READ(scratch2) rhoabso + ENDIF + ! now if we are running with one processor, also write the wave functions + ! onto that file; otherwise open new individualized file + ! + IF(ifabsorbital==1) THEN + IF(mpi_nprocs>1) THEN + IF(mpi_myproc/=0) CLOSE(UNIT=scratch2,STATUS='KEEP') + WRITE(rsfp,'(I3.3,''.'',A)') mpi_myproc,wffile(1:lenstr)//'.a' + OPEN(UNIT=scratch2,FILE=rsfp,STATUS='OLD', & + FORM='UNFORMATTED') + ENDIF + ! + DO nst = 1,nstloc + READ(scratch2) rhoescmaskOrb(:,:,:,nst) + ENDDO + END IF + CLOSE(UNIT=scratch2,STATUS='KEEP') + WRITE(6,*) 'Absorption data read from restart file' + END SUBROUTINE readabso + + PURE REAL(db) FUNCTION abs2(c) + COMPLEX(db),INTENT(IN) :: c + abs2 = REAL(c)**2 + AIMAG(c)**2 + RETURN + END FUNCTION abs2 + + +END MODULE abso_bc From 4640ac83439696ecbc4c24f373b6c728aca1b278 Mon Sep 17 00:00:00 2001 From: Kyle Godbey Date: Wed, 11 Mar 2026 15:52:34 -0400 Subject: [PATCH 3/8] Update Code/Makefile Deps Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- Code/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Makefile b/Code/Makefile index 56a8727..7a6eb10 100644 --- a/Code/Makefile +++ b/Code/Makefile @@ -126,7 +126,7 @@ dynamic.o : dynamic.f90 params.o densities.o grids.o trivial.o \ levels.o moment.o energies.o inout.o meanfield.o $(PARALLEL_SKY) \ twobody.o external.o abso_bc.o Makefile -abso_bc.o : abso_bc.f90 params.o grids.o levels.o densities.o Makefile +abso_bc.o : abso_bc.f90 params.o grids.o levels.o densities.o meanfield.o $(PARALLEL_SKY) Makefile main3d.o : main3d.f90 params.o fourier.o forces.o densities.o meanfield.o \ levels.o grids.o fragments.o $(PARALLEL_SKY) dynamic.o static.o coulomb.o user.o \ From 97a8bb4979c39d7b910125c371620b16cf3a8ce2 Mon Sep 17 00:00:00 2001 From: Kyle Godbey Date: Wed, 11 Mar 2026 15:55:47 -0400 Subject: [PATCH 4/8] Replace scalapack Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- Code/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Makefile b/Code/Makefile index 7a6eb10..ed8c06c 100644 --- a/Code/Makefile +++ b/Code/Makefile @@ -11,7 +11,7 @@ ifort ifort_seq ifort_omp : export COMPILER_SKY = ifort # # ----- define libraries----------------------------------------------------- # -seq omp debug seq_debug omp_debug ifort ifort_seq ifort_omp : export LIBS_SKY = -lfftw3 -lscalapack -lopenblas +seq omp debug seq_debug omp_debug ifort ifort_seq ifort_omp : export LIBS_SKY = -lfftw3 -llapack -lopenblas mpi mpi_debug mpi-omp : export LIBS_SKY = -lfftw3 -lscalapack -llapack -lopenblas apple apple_omp : export LIBS_SKY = -lfftw3 -framework Accelerate -L/opt/homebrew/lib # From e02d674d3bd2ed18421c479c715cc092a682b9c2 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 11 Mar 2026 19:58:36 +0000 Subject: [PATCH 5/8] Initial plan From 55a0abc558c60665230a0cfe1ef4e22d60ea8825 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 11 Mar 2026 20:02:29 +0000 Subject: [PATCH 6/8] Move absbc to per-iteration loop with density recomputation Co-authored-by: kylegodbey <5806808+kylegodbey@users.noreply.github.com> --- Code/dynamic.f90 | 18 +++++++++++++++++- Code/params.mod | Bin 0 -> 1035 bytes 2 files changed, 17 insertions(+), 1 deletion(-) create mode 100644 Code/params.mod diff --git a/Code/dynamic.f90 b/Code/dynamic.f90 index e40306b..ce2b241 100644 --- a/Code/dynamic.f90 +++ b/Code/dynamic.f90 @@ -224,7 +224,6 @@ SUBROUTINE dynamichf !$OMP END PARALLEL DO IF(tmpi) CALL collect_densities ! calculate mean fields and external fields - IF(nabsorb > 0) CALL absbc(nabsorb,iter,nt,time) CALL skyrme IF(text_timedep) CALL extfld(0.D0) CALL tinfo @@ -282,6 +281,23 @@ SUBROUTINE dynamichf !$OMP END PARALLEL DO ! sum up over nodes IF(tmpi) CALL collect_densities + ! Apply absorbing boundary conditions and recompute densities + IF(nabsorb > 0) THEN + CALL absbc(nabsorb,iter,nt,time) + rho=0.0D0 + tau=0.0D0 + current=0.0D0 + sdens=0.0D0 + sodens=0.0D0 + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst) SCHEDULE(STATIC) & + !$OMP REDUCTION(+:rho,tau,current,sdens,sodens) + DO nst=1,nstloc + CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), & + psi(:,:,:,:,nst),rho,tau,current,sdens,sodens) + ENDDO + !$OMP END PARALLEL DO + IF(tmpi) CALL collect_densities + END IF ! Step 4: eliminate center-of-mass motion if desired IF(mrescm/=0) THEN IF(MOD(iter,mrescm)==0) THEN diff --git a/Code/params.mod b/Code/params.mod new file mode 100644 index 0000000000000000000000000000000000000000..577835ac9b4a5228e8ab45f169d99c87ab261802 GIT binary patch literal 1035 zcmV+m1oZnKiwFP!000001Kn9obDKyIz0a?h8}cEk;ZxErtW3PBwUG-;*~&?UVTh>} zphb9{{QG@?_y__;W~~iT$u3XhDba7b`@NZ-n~w`?Y3U{nBKLX2=|4Q)`B6v(Q5C2Y zb35TK_2MX?TRXOc-Ft7!$Y^v$-{OZ9gDZkz%dbQQCko+{=Wn<`h1>@W1S-r|?Lqb_?5RY8O2#`VfN$bbB?=H|nE z_I(A^sqW*-oGs@Ilj@eGe_k(ESfTZ`ac`K{hB^D3J4H?-M>yR;0VUmS{qQoBN5;Iw ziJmv-soc}~vCbtj?h%&Qg&ous8uHfmR5EkV2t-Ce+(z2 zI&s6C9-LaSJgf>jEHOWKVI Date: Wed, 11 Mar 2026 20:02:50 +0000 Subject: [PATCH 7/8] Remove params.mod build artifact and add .gitignore for Code directory Co-authored-by: kylegodbey <5806808+kylegodbey@users.noreply.github.com> --- Code/.gitignore | 1 + Code/params.mod | Bin 1035 -> 0 bytes 2 files changed, 1 insertion(+) create mode 100644 Code/.gitignore delete mode 100644 Code/params.mod diff --git a/Code/.gitignore b/Code/.gitignore new file mode 100644 index 0000000..63a7748 --- /dev/null +++ b/Code/.gitignore @@ -0,0 +1 @@ +*.mod diff --git a/Code/params.mod b/Code/params.mod deleted file mode 100644 index 577835ac9b4a5228e8ab45f169d99c87ab261802..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1035 zcmV+m1oZnKiwFP!000001Kn9obDKyIz0a?h8}cEk;ZxErtW3PBwUG-;*~&?UVTh>} zphb9{{QG@?_y__;W~~iT$u3XhDba7b`@NZ-n~w`?Y3U{nBKLX2=|4Q)`B6v(Q5C2Y zb35TK_2MX?TRXOc-Ft7!$Y^v$-{OZ9gDZkz%dbQQCko+{=Wn<`h1>@W1S-r|?Lqb_?5RY8O2#`VfN$bbB?=H|nE z_I(A^sqW*-oGs@Ilj@eGe_k(ESfTZ`ac`K{hB^D3J4H?-M>yR;0VUmS{qQoBN5;Iw ziJmv-soc}~vCbtj?h%&Qg&ous8uHfmR5EkV2t-Ce+(z2 zI&s6C9-LaSJgf>jEHOWKVI Date: Wed, 11 Mar 2026 20:21:54 +0000 Subject: [PATCH 8/8] Avoid double density computation when nabsorb > 0 Co-authored-by: kylegodbey <5806808+kylegodbey@users.noreply.github.com> --- Code/dynamic.f90 | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/Code/dynamic.f90 b/Code/dynamic.f90 index ce2b241..1c714c9 100644 --- a/Code/dynamic.f90 +++ b/Code/dynamic.f90 @@ -268,27 +268,22 @@ SUBROUTINE dynamichf current=0.0D0 sdens=0.0D0 sodens=0.0D0 - ! propagate to end of step, accumulate densities + ! propagate to end of step; accumulate densities only if no absorbing boundaries !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,ps4) SCHEDULE(STATIC) & !$OMP REDUCTION(+:rho,tau,current,sdens,sodens) DO nst=1,nstloc ps4=psi(:,:,:,:,nst) CALL tstep(isospin(globalindex(nst)),mxpact,ps4) - CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), & + IF(nabsorb == 0) CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), & ps4,rho,tau,current,sdens,sodens) psi(:,:,:,:,nst)=ps4 ENDDO !$OMP END PARALLEL DO ! sum up over nodes - IF(tmpi) CALL collect_densities - ! Apply absorbing boundary conditions and recompute densities + IF(nabsorb == 0 .AND. tmpi) CALL collect_densities + ! Apply absorbing boundary conditions and compute densities from masked wavefunctions IF(nabsorb > 0) THEN CALL absbc(nabsorb,iter,nt,time) - rho=0.0D0 - tau=0.0D0 - current=0.0D0 - sdens=0.0D0 - sodens=0.0D0 !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst) SCHEDULE(STATIC) & !$OMP REDUCTION(+:rho,tau,current,sdens,sodens) DO nst=1,nstloc