From 6690843a8ed9227936202fbc0c8e6832462e789e Mon Sep 17 00:00:00 2001 From: yoricklassmann Date: Sun, 12 Jul 2026 14:14:51 +0200 Subject: [PATCH 1/3] Uncoupled SpawnModule from SelectionModule and moved FMS_RemoveDead to SelectionModule --- src/dynamics.f90 | 4 +- src/fortran_interfaces.inc | 5 - src/modules/Makefile | 4 +- src/modules/SelectionModule.f90 | 333 ++++++++++++++++++++++---------- src/modules/SpawnModule.f90 | 101 +++++++++- src/openfms.F90 | 1 + src/remove_dead.f90 | 224 --------------------- 7 files changed, 336 insertions(+), 336 deletions(-) delete mode 100644 src/remove_dead.f90 diff --git a/src/dynamics.f90 b/src/dynamics.f90 index be7e1e7..eb08a58 100644 --- a/src/dynamics.f90 +++ b/src/dynamics.f90 @@ -15,7 +15,7 @@ !! @ingroup propagation !< recursive subroutine FMS_Dynamics(Bundle, Timestep, GoalTime) - use FMSModule, only: FMS_Shutdown, FMS_RemoveDead + use FMSModule, only: FMS_Shutdown use GlobalModule, only: DefReal, DefInt, FPZero, fmiOut, & & FMS_StepRejected, FMS_PrintMessg, FMSWorkingDir, FMS_RejectStep, & & fmzWriteEveryStep, gldCoupTimeStep, & @@ -23,7 +23,7 @@ recursive subroutine FMS_Dynamics(Bundle, Timestep, GoalTime) & glzStochastic, glnStepsRejected, FMS_DieError use BundleModule use BundleIOModule, only: FMS_Output - use SelectionModule, only: FMS_StochasticCollapse + use SelectionModule, only: FMS_StochasticCollapse, FMS_RemoveDead use RestartModule, only: PutRestart use PropagationModule, only: FMS_Monitor implicit none diff --git a/src/fortran_interfaces.inc b/src/fortran_interfaces.inc index a280454..b6595c2 100644 --- a/src/fortran_interfaces.inc +++ b/src/fortran_interfaces.inc @@ -37,11 +37,6 @@ real (kind=DefReal) :: SimulationTIme end subroutine FMS_ReadNameList - subroutine FMS_RemoveDead(B1) - use BundleModule - type(T_TrajectoryBundle), intent(inout) :: B1 - end subroutine FMS_RemoveDead - subroutine FMS_Shutdown(B1, T1, terminate) use TrajectoryModule use BundleModule diff --git a/src/modules/Makefile b/src/modules/Makefile index 7c60c4d..24405e5 100644 --- a/src/modules/Makefile +++ b/src/modules/Makefile @@ -101,9 +101,9 @@ BundleIOModule.o: BundleIOModule.f90 BundleCalcsModule.o BundleModule.o \ SMDModule.o OverlapModule.o SpawnModule.o SelectionModule.o: SelectionModule.f90 GlobalModule.o RandomModule.o OverlapModule.o \ - TrajectoryModule.o BundleModule.o BundleCalcsModule.o + TrajectoryModule.o BundleModule.o BundleCalcsModule.o SpawnModule.o -SpawnModule.o: SpawnModule.f90 GlobalModule.o SelectionModule.o FMSModule.o \ +SpawnModule.o: SpawnModule.f90 GlobalModule.o FMSModule.o \ TrajectoryModule.o TrajectoryCalcsModule.o TrajectoryIOModule.o \ BundleModule.o BundleCalcsModule.o OverlapModule.o VerletModule.o diff --git a/src/modules/SelectionModule.f90 b/src/modules/SelectionModule.f90 index cb87f03..9fb2474 100644 --- a/src/modules/SelectionModule.f90 +++ b/src/modules/SelectionModule.f90 @@ -3,14 +3,15 @@ module SelectionModule use GlobalModule use BundleModule use TrajectoryModule - use BundleCalcsModule, only: FMS_bH, FMS_Norm + use BundleCalcsModule, only: FMS_bH, FMS_Norm, FMS_Mulliken use OverlapModule, only: overlap use RandomModule, only: fms_ranb + use SpawnModule, only: spdCFthresh, spawn_couple, spdpopToSpawn implicit none private public :: print_stochastic_selection_params - public :: FMS_StochasticCollapse, FMS_CalculateSelectionTime + public :: FMS_StochasticCollapse, FMS_RemoveDead ! The standard overlap threshold is 1/e, and allows us to differentiate ! betweeen normal and premature selections, where the latter are defined @@ -457,106 +458,6 @@ subroutine perform_stochastic_selection(B1, selectionTime) end subroutine perform_stochastic_selection -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine FMS_CalculateSelectionTime(parent_s, child_s, child_i) -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! -! Specific to AIMSWISS. -! Called by SpawnModule:FMS_Spawn subroutine -! -! parent_s and child_s are the parent and child TBFs at the spawning time, -! whose forces are used to calculate the selection time. Their state is not -! changed in this subroutine! -! -! child_i is the child TBF at the entry time, which will be added to the -! Bundle object after this subroutine is finished. Its state is modified -! in this subroutine, by setting its T1%SWISS%SelectionTime variable to -! its current time + the decoherence time with respect to its parent. -! -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - type(T_Trajectory), intent(in) :: parent_s, child_s - type(T_Trajectory), intent(inout) :: child_i - real(kind=DefReal) :: decoherenceTime - - decoherenceTime = FMS_CalculateDecoherenceTime(parent_s, child_s, parent_s%NumParticles) - child_i%SWISS%SelectionTime = child_i%SWISS%BirthDate + decoherenceTime - child_i%SWISS%ParentOverlap = abs(overlap(parent_s, child_s))**2 - - write (fmiOut, '(a,i0,a,f0.2)') 'SWISS: Trajectory ', parent_s%TrajID, ' and '// & - 'its child will decohere at t = ', child_i%SWISS%SelectionTime - write (fmiOut, '(a,f5.3)') 'Their current absolute overlap is ', child_i%SWISS%ParentOverlap - - end subroutine FMS_CalculateSelectionTime - -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - function FMS_CalculateDecoherenceTime(parent, child, npart) result(decoherenceTime) -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! -! Specific to AIMSWISS. -! Called by FMS_CalculateSelectionTime: -! -! Caclculate the decoherence time (tau_D) in the following way: -! -! First, we calculate the decoherence rate (Gamma_D) rate via -! -! Gamma_D = (F_P - F_C)^T * alpha^-1 * (F_P - F_C) / 4, -! -! where F_P and F_C are the gradients acting on the parent and child TBF, -! respectively, and alpha is a matrix containing the widths of the TBFs. -! -! Second, we determine tau_D by taking the square root of Gamma_D and -! its reciprocal -! -! tau_D = 1. / sqrt(Gamma_D) -! -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - type(T_Trajectory), intent(in) :: parent, child - integer(kind=DefInt), intent(in) :: npart - real(kind=DefReal) :: decoherenceTime - integer(kind=DefInt) :: ipart, jdim, ishifted, istart, ndim, & - parentSt, childSt - real(kind=DefReal) :: forceElementDiff2, decoherenceRate2, & - parentForceElement, childForceElement, & - summand2, width - -! Calculation of the squared decoherence rate - parentSt = parent%StateID - childSt = child%StateID -! stochastic selection is currently only implemented for -! multi-state spawning (may change in the future) - if (childSt == parentSt) then - write (fmiOut, *) "SWISS: States don't differ, decoherence not possible" - call FMS_DieError('ERROR in FMS_CalculateDecoherenceTime') - end if - - decoherenceRate2 = 0.d0 - forceElementDiff2 = 0.d0 - do ipart = 1, npart - width = parent%Particle(ipart)%Width - ndim = parent%Particle(ipart)%NumDimensions - - istart = (ipart - 1) * ndim - - do jdim = 1, ndim - ishifted = jdim + istart - parentForceElement = parent%ElecStruc%DerivMat(parentSt, & - parentSt, ishifted) - childForceElement = child%ElecStruc%DerivMat(childSt, & - childSt, ishifted) - forceElementDiff2 = (parentForceElement - & - childForceElement)**2 - summand2 = forceElementDiff2 / (4.d0 * width) - decoherenceRate2 = decoherenceRate2 + summand2 - end do - end do - - decoherenceTime = 1.d0 / sqrt(decoherenceRate2) - - return - - end function FMS_CalculateDecoherenceTime - ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine FMS_CheckSelectionTime(B1, performSelection, selectionTime) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -1410,4 +1311,232 @@ subroutine copy_state_bundles_to_original_bundle(B1, BundleSS) end subroutine copy_state_bundles_to_original_bundle +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine FMS_BuildCoupled_RemoveDead(B1, Coupled) +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + type(T_TrajectoryBundle), intent(in) :: B1 + integer(kind=DefInt), dimension(B1%NumTraj, B1%NumTraj), & + intent(inout) :: Coupled + integer(kind=DefInt) :: ntraj, i, j + + ntraj = B1%NumTraj + + do i = 2, ntraj + do j = 1, i + if (abs(FMS_bH(B1, i, j)) > FPZero) then + Coupled(i, j) = 1 + Coupled(j, i) = 1 + end if + end do + end do + +! Uncover all indirect connections between TBFs + call FMS_ConvergeCoupled(ntraj, Coupled) + + end subroutine FMS_BuildCoupled_RemoveDead + + subroutine FMS_Count_Nr_Dead(B1, Coupled, nr_dead) + type(T_TrajectoryBundle), intent(inout) :: B1 + integer(kind=DefInt), dimension(B1%NumTraj, B1%NumTraj), & + intent(inout) :: Coupled + integer(kind=DefInt), intent(out) :: nr_dead + + real(DefReal) :: Population(B1%NumTraj) + logical :: NotCoupled, OnIgnoreState, PopBelowThresh, & + MarkForDeath, OverSpawnThresh, ForceDead + ! GAIMS added + real(kind=DefReal) :: pop + ! GAIMS added end + + integer(DefInt) :: ntraj, i, TrajID, StateID, nstate, n + + ntraj = B1%NumTraj + nstate = B1%NumStates + + Population = abs(FMS_Mulliken(B1)) + nr_dead = 0 + do i = 1, ntraj + + TrajID = B1%Trajectory(i)%TrajID + StateID = B1%Trajectory(i)%StateID + + OverSpawnThresh = .false. + do n = 1, nstate + OverSpawnThresh = OverSpawnThresh .or. spawn_couple(B1%Trajectory(i), n) > spdCFThresh + end do + + ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + ! Update DeadTime for those who will not be killed + NotCoupled = all(Coupled(:, i) == 0) + OnIgnoreState = (StateID == glIgnoreState) + PopBelowThresh = (Population(i) < spdPopToSpawn) + ! GAIMS changed + pop = 0.d0 + do n = 1, ntraj + if (B1%Trajectory(n)%CBF == B1%Trajectory(i)%CBF) then + pop = pop + Population(n) + end if + end do + PopBelowThresh = (pop < spdPopToSpawn) + ! GAIMS changed end + MarkForDeath = ((OnIgnoreState .and. NotCoupled .and. .not. OverSpawnThresh) & + .or. (PopBelowThresh .and. NotCoupled)) + ForceDead = any(gliForceKill(:) == TrajID) + if (ForceDead) then + B1%Trajectory(i)%DeadTime = -9999.0d0 + end if + MarkForDeath = (MarkForDeath .or. ForceDead) + + if (.not. MarkForDeath) then + B1%Trajectory(i)%DeadTime = B1%CurrentTime + else + write (fmiOut, '(a,i0,a)') 'Traj ', TrajID, ' was marked for death' + end if + + ! Workout who is getting killed + if (B1%Trajectory(i)%is_dead()) then + nr_dead = nr_dead + 1 + if (ForceDead) then + write (fmiOut, '(a,i0,a,i0)') '** Force Killing trajectory ', TrajID, ' on state ', StateID + else if (OnIgnoreState) then + write (fmiOut, '(a,i0,a,i0)') '** Killing trajectory ', TrajID, ' on state ', StateID + else + write (fmiOut, '(a,i0,a,f6.5)') '** Killing trajectory ', TrajID, ' pop < ', spdPopToSpawn + end if + end if + end do + + end subroutine FMS_Count_Nr_Dead + +! Copyright Todd J. Martinez and Raphael D. Levine, 1994 +!> +!! Mark any trajectories to be removed as 'dead', then remove them +!! and rescale the bundle +!! +!! Trajectories are marked for deletion if they are on IgnoreState, +!! or have population less than PopToSpawn. If either of these is +!! true, then a CountDown timer is started for that trajectory. If, +!! after 5 fs, one of those conditions is still true, the trajectory +!! is removed from the Bundle%Trajectory array and appended to the +!! Bundle%DeadTraj array. The bundle matrices are rescaled +!! accordingly. +!! +!! @ingroup propagation +!< +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine FMS_RemoveDead(B1) + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + type(T_TrajectoryBundle), intent(inout) :: B1 + type(T_TrajectoryBundle) :: BTemp + + integer(kind=DefInt) :: iTraj, nDead, iLive, iDead, jTraj + integer(kind=DefInt) :: iCent, jCent, jDead, jLive + integer(kind=DefInt) :: iCBF, jCBF, nCBF, n2CBF + + integer(DefInt), dimension(B1%NumTraj, B1%NumTraj) :: Coupled + + integer(DefInt) :: ntraj, i + + ntraj = B1%NumTraj + + Coupled = 0 + + call FMS_BuildCoupled_RemoveDead(B1, Coupled) + + do i = 1, ntraj + Coupled(i, i) = 0 + end do + + call FMS_Count_Nr_Dead(B1, Coupled, ndead) + + ! Remove dead trajectories + if (ndead > 0) then + call BTemp%create(numtraj=B1%NumTraj - nDead, & + numdeadtraj=B1%NumDeadTraj + nDead, & + numstates=B1%NumStates, & + numparticles=B1%NumParticles, & + ncbfs=B1%NCBFs) + BTemp%CurrentTime = B1%CurrentTime + + ! Copy old dead trajectories to new bundle + do iTraj = 1, B1%NumDeadTraj + call BTemp%DeadTraj(iTraj)%copy_from(B1%DeadTraj(iTraj)) + do jTraj = 1, B1%NumDeadTraj + BTemp%DeadH(iTraj, jTraj) = B1%DeadH(iTraj, jTraj) + end do + end do + + iLive = 0 + iDead = 0 + nCBF = 0 + write (fmiOut, '(a,/,a)') 'Reconstructing trajectory bundle. Living trajectories:' + + do iTraj = 1, B1%NumTraj + if (.not. B1%Trajectory(iTraj)%is_dead()) then + ! Alive: copy trajectory to new bundle + iLive = iLive + 1 + call BTemp%Trajectory(iLive)%copy_from(B1%Trajectory(iTraj)) + ! Get the CBF identifier. Here we use the fact that + ! the 3 triplet trajecories are always stored in + ! successive order starting with Ms=2 + if (BTemp%Trajectory(iLive)%Ms == 2) then + nCBF = nCBF + 1 + end if + if (B1%Trajectory(iTraj)%triplet) then + write (fmiOut, '(4X,I4," S=1 Ms=",I0," CBF ",I0)') iTraj, BTemp%Trajectory(iLive)%Ms - 2, nCBF + else + write (fmiOut, '(4X,I4," S=0 Ms=",I0," CBF ",I0)') iTraj, 0, nCBF + end if + BTemp%Trajectory(iLive)%CBF = nCBF + ! Copy centroids over too + ![bfec + if (glzCentroids) then + jLive = 0 + if (BTemp%Trajectory(iLive)%Ms == 2) then + n2CBF = 0 + do jTraj = 1, iTraj - 1 + if (.not. B1%Trajectory(jTraj)%is_dead() .and. (B1%Trajectory(jTraj)%Ms == 2)) then + n2CBF = n2CBF + 1 + iCBF = B1%Trajectory(iTraj)%CBF + jCBF = B1%Trajectory(jTraj)%CBF + iCent = ((iCBF - 2) * (iCBF - 1)) / 2 + jCBF + jCent = ((nCBF - 2) * (nCBF - 1)) / 2 + n2CBF + call BTemp%Centroids(jCent)%copy_from(B1%Centroids(iCent)) + BTemp%Centroids(jCent)%CentID(1) = nCBF + BTemp%Centroids(jCent)%CentID(2) = n2CBF + end if + end do + end if + end if + + else + ! Dead: add the trajectory to the graveyard + iDead = iDead + 1 + ! TODO: Use FMS_AssignTrajectory explicitly + ! call FMS_AssignTrajectory(BTemp%DeadTraj(B1%NumDeadTraj + iDead), B1%Trajectory(iTraj)) + call BTemp%DeadTraj(B1%NumDeadTraj + iDead)%copy_from(B1%Trajectory(iTraj)) + + ! Set dead time of recently killed trajectory to be current time + ! to make it easier for user to restart killed trajectories! + BTemp%DeadTraj(B1%NumDeadTraj + iDead)%DeadTime = BTemp%CurrentTime + ! Copy new hamiltonian elements into dead hamiltonian + jDead = 0 + do jTraj = 1, iTraj + if (B1%Trajectory(jTraj)%is_dead()) then + jDead = jDead + 1 + BTemp%DeadH(B1%NumDeadTraj + iDead, B1%NumDeadTraj + jDead) = FMS_bH(B1, iTraj, jTraj) + BTemp%DeadH(B1%NumDeadTraj + jDead, B1%NumDeadTraj + iDead) = FMS_bH(B1, jTraj, iTraj) + end if + end do + end if + end do + + BTemp%NCBFs = nCBF + call B1%copy_from(BTemp) + call BTemp%destroy() + + end if + + end subroutine FMS_RemoveDead + end module SelectionModule diff --git a/src/modules/SpawnModule.f90 b/src/modules/SpawnModule.f90 index bd50128..ea18f35 100644 --- a/src/modules/SpawnModule.f90 +++ b/src/modules/SpawnModule.f90 @@ -102,7 +102,6 @@ subroutine FMS_Spawn(B1, TimeStep) ! child_i <--- child_s -----> child_f ! use ElecStrucModule - use SelectionModule, only: FMS_CalculateSelectionTime type(T_TrajectoryBundle), target, intent(inout) :: B1 real(kind=DefReal), intent(in) :: TimeStep @@ -1330,4 +1329,104 @@ function FMS_AdjustEnergy(TChild, TParent, ScaleVector) result(Success) end function FMS_AdjustEnergy +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine FMS_CalculateSelectionTime(parent_s, child_s, child_i) +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! +! Specific to AIMSWISS. +! Called by SpawnModule:FMS_Spawn subroutine +! +! parent_s and child_s are the parent and child TBFs at the spawning time, +! whose forces are used to calculate the selection time. Their state is not +! changed in this subroutine! +! +! child_i is the child TBF at the entry time, which will be added to the +! Bundle object after this subroutine is finished. Its state is modified +! in this subroutine, by setting its T1%SWISS%SelectionTime variable to +! its current time + the decoherence time with respect to its parent. +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + type(T_Trajectory), intent(in) :: parent_s, child_s + type(T_Trajectory), intent(inout) :: child_i + real(kind=DefReal) :: decoherenceTime + + decoherenceTime = FMS_CalculateDecoherenceTime(parent_s, child_s, parent_s%NumParticles) + child_i%SWISS%SelectionTime = child_i%SWISS%BirthDate + decoherenceTime + child_i%SWISS%ParentOverlap = abs(overlap(parent_s, child_s))**2 + + write (fmiOut, '(a,i0,a,f0.2)') 'SWISS: Trajectory ', parent_s%TrajID, ' and '// & + 'its child will decohere at t = ', child_i%SWISS%SelectionTime + write (fmiOut, '(a,f5.3)') 'Their current absolute overlap is ', child_i%SWISS%ParentOverlap + + end subroutine FMS_CalculateSelectionTime + +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + function FMS_CalculateDecoherenceTime(parent, child, npart) result(decoherenceTime) +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! +! Specific to AIMSWISS. +! Called by FMS_CalculateSelectionTime: +! +! Caclculate the decoherence time (tau_D) in the following way: +! +! First, we calculate the decoherence rate (Gamma_D) rate via +! +! Gamma_D = (F_P - F_C)^T * alpha^-1 * (F_P - F_C) / 4, +! +! where F_P and F_C are the gradients acting on the parent and child TBF, +! respectively, and alpha is a matrix containing the widths of the TBFs. +! +! Second, we determine tau_D by taking the square root of Gamma_D and +! its reciprocal +! +! tau_D = 1. / sqrt(Gamma_D) +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + type(T_Trajectory), intent(in) :: parent, child + integer(kind=DefInt), intent(in) :: npart + real(kind=DefReal) :: decoherenceTime + integer(kind=DefInt) :: ipart, jdim, ishifted, istart, ndim, & + parentSt, childSt + real(kind=DefReal) :: forceElementDiff2, decoherenceRate2, & + parentForceElement, childForceElement, & + summand2, width + +! Calculation of the squared decoherence rate + parentSt = parent%StateID + childSt = child%StateID +! stochastic selection is currently only implemented for +! multi-state spawning (may change in the future) + if (childSt == parentSt) then + write (fmiOut, *) "SWISS: States don't differ, decoherence not possible" + call FMS_DieError('ERROR in FMS_CalculateDecoherenceTime') + end if + + decoherenceRate2 = 0.d0 + forceElementDiff2 = 0.d0 + do ipart = 1, npart + width = parent%Particle(ipart)%Width + ndim = parent%Particle(ipart)%NumDimensions + + istart = (ipart - 1) * ndim + + do jdim = 1, ndim + ishifted = jdim + istart + parentForceElement = parent%ElecStruc%DerivMat(parentSt, & + parentSt, ishifted) + childForceElement = child%ElecStruc%DerivMat(childSt, & + childSt, ishifted) + forceElementDiff2 = (parentForceElement - & + childForceElement)**2 + summand2 = forceElementDiff2 / (4.d0 * width) + decoherenceRate2 = decoherenceRate2 + summand2 + end do + end do + + decoherenceTime = 1.d0 / sqrt(decoherenceRate2) + + return + + end function FMS_CalculateDecoherenceTime + end module SpawnModule diff --git a/src/openfms.F90 b/src/openfms.F90 index afd11c4..136e042 100644 --- a/src/openfms.F90 +++ b/src/openfms.F90 @@ -21,6 +21,7 @@ program OpenFMS use ElecStrucModule, only: FMS_ESInit use RestartModule, only: inIRestart, RestartTime, getRestart use PropagationModule, only: FMS_SetTimeStep + use SelectionModule, only: FMS_RemoveDead implicit none type(T_TrajectoryBundle) :: Bundle diff --git a/src/remove_dead.f90 b/src/remove_dead.f90 deleted file mode 100644 index d26af5a..0000000 --- a/src/remove_dead.f90 +++ /dev/null @@ -1,224 +0,0 @@ -! Copyright Todd J. Martinez and Raphael D. Levine, 1994 -!> -!! Mark any trajectories to be removed as 'dead', then remove them -!! and rescale the bundle -!! -!! Trajectories are marked for deletion if they are on IgnoreState, -!! or have population less than PopToSpawn. If either of these is -!! true, then a CountDown timer is started for that trajectory. If, -!! after 5 fs, one of those conditions is still true, the trajectory -!! is removed from the Bundle%Trajectory array and appended to the -!! Bundle%DeadTraj array. The bundle matrices are rescaled -!! accordingly. -!! -!! @ingroup propagation -!< -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -subroutine FMS_RemoveDead(B1) - ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - use GlobalModule, only: defInt, defReal, fmiOut, FPZero, & - glIgnoreState, glzCentroids, gliForceKill - use BundleModule - use BundleCalcsModule, only: FMS_bH, FMS_Mulliken - use SpawnModule, only: spdCFthresh, spawn_couple, spdpopToSpawn - implicit none - type(T_TrajectoryBundle), intent(inout) :: B1 - type(T_TrajectoryBundle) :: BTemp - - integer, parameter :: MaxIter = 50 - integer(kind=DefInt) :: iTraj, nDead, iLive, iDead, jTraj - integer(kind=DefInt) :: iCent, jCent, jDead, jLive - integer(kind=DefInt) :: iCBF, jCBF, nCBF, n2CBF - - integer(DefInt), dimension(B1%NumTraj, B1%NumTraj) :: Coupled, Coupled_prev - - real(DefReal) :: Population(B1%NumTraj) - -! GAIMS added - real(kind=DefReal) :: pop -! GAIMS added end - - integer(DefInt) :: ntraj, i, j, TrajID, StateID, nstate, n - - logical :: NotCoupled, OnIgnoreState, PopBelowThresh, & - MarkForDeath, OverSpawnThresh, ForceDead - - ntraj = B1%NumTraj - nstate = B1%NumStates - -! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . -! First, decompose the hamiltonian into a block diagonal representation - ! work out the coupling matrix - Coupled = 0 - do i = 2, ntraj - do j = 1, i - if (abs(FMS_bH(B1, i, j)) > FPZero) then - Coupled(i, j) = 1 - Coupled(j, i) = 1 - end if - end do - end do - -! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . -! This matrix has the properties that it is non-zero is -! Coup^1 : directy connected -! Coup^2 : connected by 1 or less common trajectory -! Coup^3 : connected by 2 or less common trajectories -! We will iterate the matrix multiplication to convergence - do - Coupled_prev = Coupled - - Coupled = min(matmul(Coupled, Coupled), 1) - - if (all(Coupled == Coupled_prev)) exit - end do - -! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . -! Replace the diagonal with zero - do i = 1, ntraj - Coupled(i, i) = 0 - end do - - Population = abs(FMS_Mulliken(B1)) - - ndead = 0 - do i = 1, ntraj - - TrajID = B1%Trajectory(i)%TrajID - StateID = B1%Trajectory(i)%StateID - - OverSpawnThresh = .false. - do n = 1, nstate - OverSpawnThresh = OverSpawnThresh .or. spawn_couple(B1%Trajectory(i), n) > spdCFThresh - end do - -! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . -! Update DeadTime for those who will not be killed - NotCoupled = all(Coupled(:, i) == 0) - OnIgnoreState = (StateID == glIgnoreState) - PopBelowThresh = (Population(i) < spdPopToSpawn) -! GAIMS changed - pop = 0.d0 - do n = 1, ntraj - if (B1%Trajectory(n)%CBF == B1%Trajectory(i)%CBF) then - pop = pop + Population(n) - end if - end do - PopBelowThresh = (pop < spdPopToSpawn) -! GAIMS changed end - MarkForDeath = ((OnIgnoreState .and. NotCoupled .and. .not. OverSpawnThresh) & - .or. (PopBelowThresh .and. NotCoupled)) - ForceDead = any(gliForceKill(:) == TrajID) - if (ForceDead) then - B1%Trajectory(i)%DeadTime = -9999.0d0 - end if - MarkForDeath = (MarkForDeath .or. ForceDead) - - if (.not. MarkForDeath) then - B1%Trajectory(i)%DeadTime = B1%CurrentTime - else - write (fmiOut, '(a,i0,a)') 'Traj ', TrajID, ' was marked for death' - end if - -! Workout who is getting killed - if (B1%Trajectory(i)%is_dead()) then - ndead = ndead + 1 - if (ForceDead) then - write (fmiOut, '(a,i0,a,i0)') '** Force Killing trajectory ', TrajID, ' on state ', StateID - else if (OnIgnoreState) then - write (fmiOut, '(a,i0,a,i0)') '** Killing trajectory ', TrajID, ' on state ', StateID - else - write (fmiOut, '(a,i0,a,f6.5)') '** Killing trajectory ', TrajID, ' pop < ', spdPopToSpawn - end if - end if - end do - -! Remove dead trajectories - if (ndead > 0) then - call BTemp%create(numtraj=B1%NumTraj - nDead, & - numdeadtraj=B1%NumDeadTraj + nDead, & - numstates=B1%NumStates, & - numparticles=B1%NumParticles, & - ncbfs=B1%NCBFs) - BTemp%CurrentTime = B1%CurrentTime - -! Copy old dead trajectories to new bundle - do iTraj = 1, B1%NumDeadTraj - call BTemp%DeadTraj(iTraj)%copy_from(B1%DeadTraj(iTraj)) - do jTraj = 1, B1%NumDeadTraj - BTemp%DeadH(iTraj, jTraj) = B1%DeadH(iTraj, jTraj) - end do - end do - - iLive = 0 - iDead = 0 - nCBF = 0 - write (fmiOut, '(a,/,a)') 'Reconstructing trajectory bundle. Living trajectories:' - - do iTraj = 1, B1%NumTraj - if (.not. B1%Trajectory(iTraj)%is_dead()) then -! Alive: copy trajectory to new bundle - iLive = iLive + 1 - call BTemp%Trajectory(iLive)%copy_from(B1%Trajectory(iTraj)) -! Get the CBF identifier. Here we use the fact that -! the 3 triplet trajecories are always stored in -! successive order starting with Ms=2 - if (BTemp%Trajectory(iLive)%Ms == 2) then - nCBF = nCBF + 1 - end if - if (B1%Trajectory(iTraj)%triplet) then - write (fmiOut, '(4X,I4," S=1 Ms=",I0," CBF ",I0)') iTraj, BTemp%Trajectory(iLive)%Ms - 2, nCBF - else - write (fmiOut, '(4X,I4," S=0 Ms=",I0," CBF ",I0)') iTraj, 0, nCBF - end if - BTemp%Trajectory(iLive)%CBF = nCBF -! Copy centroids over too -![bfec - if (glzCentroids) then - jLive = 0 - if (BTemp%Trajectory(iLive)%Ms == 2) then - n2CBF = 0 - do jTraj = 1, iTraj - 1 - if (.not. B1%Trajectory(jTraj)%is_dead() .and. (B1%Trajectory(jTraj)%Ms == 2)) then - n2CBF = n2CBF + 1 - iCBF = B1%Trajectory(iTraj)%CBF - jCBF = B1%Trajectory(jTraj)%CBF - iCent = ((iCBF - 2) * (iCBF - 1)) / 2 + jCBF - jCent = ((nCBF - 2) * (nCBF - 1)) / 2 + n2CBF - call BTemp%Centroids(jCent)%copy_from(B1%Centroids(iCent)) - BTemp%Centroids(jCent)%CentID(1) = nCBF - BTemp%Centroids(jCent)%CentID(2) = n2CBF - end if - end do - end if - end if - - else -! Dead: add the trajectory to the graveyard - iDead = iDead + 1 - ! TODO: Use FMS_AssignTrajectory explicitly - ! call FMS_AssignTrajectory(BTemp%DeadTraj(B1%NumDeadTraj + iDead), B1%Trajectory(iTraj)) - call BTemp%DeadTraj(B1%NumDeadTraj + iDead)%copy_from(B1%Trajectory(iTraj)) - -! Set dead time of recently killed trajectory to be current time -! to make it easier for user to restart killed trajectories! - BTemp%DeadTraj(B1%NumDeadTraj + iDead)%DeadTime = BTemp%CurrentTime -! Copy new hamiltonian elements into dead hamiltonian - jDead = 0 - do jTraj = 1, iTraj - if (B1%Trajectory(jTraj)%is_dead()) then - jDead = jDead + 1 - BTemp%DeadH(B1%NumDeadTraj + iDead, B1%NumDeadTraj + jDead) = FMS_bH(B1, iTraj, jTraj) - BTemp%DeadH(B1%NumDeadTraj + jDead, B1%NumDeadTraj + iDead) = FMS_bH(B1, jTraj, iTraj) - end if - end do - end if - end do - - BTemp%NCBFs = nCBF - call B1%copy_from(BTemp) - call BTemp%destroy() - - end if - -end subroutine FMS_RemoveDead From 6f4cb2d0aafa54bf762d3393a12b0158e207cf05 Mon Sep 17 00:00:00 2001 From: yoricklassmann Date: Sun, 12 Jul 2026 14:30:17 +0200 Subject: [PATCH 2/3] Forgot to remove remove_dead.f90 from OpenFMSSources.cmake --- cmake/OpenFMSSources.cmake | 1 - 1 file changed, 1 deletion(-) diff --git a/cmake/OpenFMSSources.cmake b/cmake/OpenFMSSources.cmake index b042a74..1cadec8 100644 --- a/cmake/OpenFMSSources.cmake +++ b/cmake/OpenFMSSources.cmake @@ -36,6 +36,5 @@ set(OPENFMS_CORE_SOURCES src/particle_types.f90 src/read_geometry.f90 src/read_namelist.f90 - src/remove_dead.f90 src/shutdown.f90 ) From 06683a454bb690a90b3df3f76900ba9fd028b8f2 Mon Sep 17 00:00:00 2001 From: yoricklassmann Date: Sun, 12 Jul 2026 14:37:49 +0200 Subject: [PATCH 3/3] Updated SelectionModule to use spawn_params derived type. --- src/modules/SelectionModule.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/modules/SelectionModule.f90 b/src/modules/SelectionModule.f90 index 9fb2474..e5eb7bc 100644 --- a/src/modules/SelectionModule.f90 +++ b/src/modules/SelectionModule.f90 @@ -6,7 +6,7 @@ module SelectionModule use BundleCalcsModule, only: FMS_bH, FMS_Norm, FMS_Mulliken use OverlapModule, only: overlap use RandomModule, only: fms_ranb - use SpawnModule, only: spdCFthresh, spawn_couple, spdpopToSpawn + use SpawnModule, only: spawn_params, spawn_couple implicit none private @@ -1362,14 +1362,14 @@ subroutine FMS_Count_Nr_Dead(B1, Coupled, nr_dead) OverSpawnThresh = .false. do n = 1, nstate - OverSpawnThresh = OverSpawnThresh .or. spawn_couple(B1%Trajectory(i), n) > spdCFThresh + OverSpawnThresh = OverSpawnThresh .or. spawn_couple(B1%Trajectory(i), n) > spawn_params%CFThresh end do ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ! Update DeadTime for those who will not be killed NotCoupled = all(Coupled(:, i) == 0) OnIgnoreState = (StateID == glIgnoreState) - PopBelowThresh = (Population(i) < spdPopToSpawn) + PopBelowThresh = (Population(i) < spawn_params%PopToSpawn) ! GAIMS changed pop = 0.d0 do n = 1, ntraj @@ -1377,7 +1377,7 @@ subroutine FMS_Count_Nr_Dead(B1, Coupled, nr_dead) pop = pop + Population(n) end if end do - PopBelowThresh = (pop < spdPopToSpawn) + PopBelowThresh = (pop < spawn_params%PopToSpawn) ! GAIMS changed end MarkForDeath = ((OnIgnoreState .and. NotCoupled .and. .not. OverSpawnThresh) & .or. (PopBelowThresh .and. NotCoupled)) @@ -1401,7 +1401,7 @@ subroutine FMS_Count_Nr_Dead(B1, Coupled, nr_dead) else if (OnIgnoreState) then write (fmiOut, '(a,i0,a,i0)') '** Killing trajectory ', TrajID, ' on state ', StateID else - write (fmiOut, '(a,i0,a,f6.5)') '** Killing trajectory ', TrajID, ' pop < ', spdPopToSpawn + write (fmiOut, '(a,i0,a,f6.5)') '** Killing trajectory ', TrajID, ' pop < ', spawn_params%PopToSpawn end if end if end do