From 83a912bb8662fb2a8a45335d3cdc0b1091d88a3b Mon Sep 17 00:00:00 2001 From: verab98 Date: Mon, 16 Mar 2026 11:29:38 +0000 Subject: [PATCH 01/10] Avoid initialising on a triplet state Trajectories starting from a triplet are currently not initialised as a multiplet. So, starting from a triplet should not be possible until this is fixed. --- src/openfms.F90 | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/openfms.F90 b/src/openfms.F90 index 70f4475..f94b3fb 100644 --- a/src/openfms.F90 +++ b/src/openfms.F90 @@ -12,7 +12,7 @@ program OpenFMS use FMSModule use GlobalModule, only: DefInt, DefReal, fmiOut, & FMS_DeleteFile, FmsWorkingDir, timei, FMS_DieError, & - fmzWriteEveryStep, gldTimeStep, glzMinSearch + fmzWriteEveryStep, gldTimeStep, glzMinSearch, NSing, NTrip use RandomModule, only: FMS_ranb, initialize_fortran_prng use BundleModule use BundleCalcsModule, only: FMS_UpdateBundle @@ -417,6 +417,15 @@ subroutine initialize_bundle(B1, NumTraj, NumParticles, NumStates, InitialState) T1(ITraj)%Amplitude = DMTemp end do +! Check if we are initialising on a triplet state + if (NTrip /= 0) then + do ITraj = 1, NumTraj + if (T1(ITraj)%StateID > NSing) then + call FMS_DieError('Trajectories cannot currently start from a triplet state') + end if + end do + end if + ! Copy trajectories into bundle and clean up do ITraj = 1, NumTraj call B1%add_traj(T1(iTraj)) From af4151eb0f7936d2a6b53237a7553eff7bfedb91 Mon Sep 17 00:00:00 2001 From: verab98 Date: Mon, 16 Mar 2026 11:44:31 +0000 Subject: [PATCH 02/10] Fixing initial Amp with more than one initial TBF Before, all TBFs would start with Amp = 1, leading to a norm > 1 Changed this by adjusting Amp initialisation, but is that the right way to go? --- src/openfms.F90 | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/src/openfms.F90 b/src/openfms.F90 index f94b3fb..48baf96 100644 --- a/src/openfms.F90 +++ b/src/openfms.F90 @@ -12,7 +12,8 @@ program OpenFMS use FMSModule use GlobalModule, only: DefInt, DefReal, fmiOut, & FMS_DeleteFile, FmsWorkingDir, timei, FMS_DieError, & - fmzWriteEveryStep, gldTimeStep, glzMinSearch, NSing, NTrip + fmzWriteEveryStep, gldTimeStep, glzMinSearch, NSing, NTrip, & + NumInitBasis use RandomModule, only: FMS_ranb, initialize_fortran_prng use BundleModule use BundleCalcsModule, only: FMS_UpdateBundle @@ -407,15 +408,25 @@ subroutine initialize_bundle(B1, NumTraj, NumParticles, NumStates, InitialState) call FMS_ParticleTypes(T1(1)) +! Initialising amplitudes depending on number of initial TBFs T1%StateID = InitialState DMTemp = 1.0 - T1(1)%Amplitude = DMTemp - do ITraj = 2, NumTraj - T1(ITraj) = T1(1) - ! DH: This seems weird??? -! T1(ITraj)%Amplitude=0.0 - T1(ITraj)%Amplitude = DMTemp - end do + if (NumInitBasis > 1) then + + T1(1)%Amplitude=DMTemp + !T1(1)%Amplitude=DMTemp / NumInitBasis ! for equally distributing initial amp + + do ITraj = 2, NumTraj + T1(ITraj) = T1(1) + ! DH: This seems weird??? + T1(ITraj)%Amplitude=0.0 +! T1(ITraj)%Amplitude=DMTemp / NumInitBasis ! for equally distributing initial amp +! T1(ITraj)%Amplitude = DMTemp + end do + + else + T1(1)%Amplitude = DMTemp + endif ! Check if we are initialising on a triplet state if (NTrip /= 0) then From cf5cec5e113027767b1a074da91b317fc96e28fa Mon Sep 17 00:00:00 2001 From: verab98 Date: Mon, 16 Mar 2026 12:01:49 +0000 Subject: [PATCH 03/10] Useful print statements for checking if Bundle with 2 initial TBFs is correctly created but nothing else --- src/openfms.F90 | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/openfms.F90 b/src/openfms.F90 index 48baf96..23e5271 100644 --- a/src/openfms.F90 +++ b/src/openfms.F90 @@ -16,7 +16,7 @@ program OpenFMS NumInitBasis use RandomModule, only: FMS_ranb, initialize_fortran_prng use BundleModule - use BundleCalcsModule, only: FMS_UpdateBundle + use BundleCalcsModule, only: FMS_UpdateBundle, FMS_Norm use BundleIOModule, only: FMS_Output use InitialModule, only: iniRndSeed, inInitState use ElecStrucModule, only: FMS_ESInit @@ -442,6 +442,11 @@ subroutine initialize_bundle(B1, NumTraj, NumParticles, NumStates, InitialState) call B1%add_traj(T1(iTraj)) end do + write (fmiout, *) "Trajectories in initially in Bundle: (after adding them)", B1%NumTraj + write (fmiout, *) "And their IDs: " + write (fmiout, *) "TrajID 1:", B1%Trajectory(1)%TrajID + write (fmiout, *) "TrajID 2:", B1%Trajectory(2)%TrajID + do IParticle = 1, NumParticles call Particle(IParticle)%destroy() end do @@ -452,6 +457,16 @@ subroutine initialize_bundle(B1, NumTraj, NumParticles, NumStates, InitialState) end do deallocate (T1) + write (fmiout, *) "Trajectories in initially in Bundle: (after clean up step)", B1%NumTraj + write (fmiout, *) "And their IDs: " + write (fmiout, *) "TrajID 1:", B1%Trajectory(1)%TrajID + write (fmiout, *) "TrajID 2:", B1%Trajectory(2)%TrajID + write (fmiout, *) "And their Amplitudes: " + write (fmiout, *) "TrajID 1:", B1%Trajectory(1)%Amplitude + write (fmiout, *) "TrajID 2:", B1%Trajectory(2)%Amplitude + write (fmiout, *) "Norm of the Bundle: " + write (fmiout, *) FMS_Norm(B1) + end subroutine initialize_bundle subroutine get_cmdline() From 396af334a36ac1e9cab8adb9a5529e13adef7e66 Mon Sep 17 00:00:00 2001 From: verab98 Date: Mon, 16 Mar 2026 17:56:51 +0000 Subject: [PATCH 04/10] More useful print statements In OverlapModule for checking bundle set up --- src/modules/OverlapModule.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/modules/OverlapModule.f90 b/src/modules/OverlapModule.f90 index 74c9a3c..77cce87 100644 --- a/src/modules/OverlapModule.f90 +++ b/src/modules/OverlapModule.f90 @@ -672,6 +672,8 @@ function overlap_V_trajectory(T_i, T_j, T_c, S_ij_precalc) result(V) ! check that the correct centroid was passed ! cent_match = (it==ic .and. jt==jc) .or. (jt==ic .and. it==jc) cent_match = (CBFi == ic .and. CBFj == jc) .or. (CBFj == ic .and. CBFi == jc) + write (fmiout, *) "(CBFij == ijc .and. CBFji == jic)", "CBFi, ic, CBFj, jc", CBFi, ic, CBFj, jc + flush (fmiout) ! cent_match = (i==ic .and. j==jc) .or. (j==ic .and. i==jc) if (.not. cent_match) then call FMS_DieError('overlap_V_trajectory :: centroids dont match') From c2a7a009c64df3dd619b50850f226126849497ed Mon Sep 17 00:00:00 2001 From: verab98 Date: Tue, 7 Apr 2026 16:35:37 +0100 Subject: [PATCH 05/10] (WiP) Separated Stochastic Selection Changes to SelectionModule: - alternative selection if singlet and triplet states - selection only happens for singlet only or triplet only blocks of coupled trajectories (not for S/T mixed blocks) Other changes: - Starting with more than one initial Gaussian (NumInitBasis > 1) - Displacing Gaussians using swarm in SamplingModule (this only works for the SOC model, as it assumes 1D) --- src/dynamics.f90 | 24 + src/modules/BundleCalcsModule.f90 | 38 +- src/modules/OverlapModule.f90 | 7 +- src/modules/SamplingModule.f90 | 86 ++-- src/modules/SelectionModule.f90 | 730 +++++++++++++++++++++++++++++- src/openfms.F90 | 43 +- 6 files changed, 877 insertions(+), 51 deletions(-) diff --git a/src/dynamics.f90 b/src/dynamics.f90 index d605669..b6a038d 100644 --- a/src/dynamics.f90 +++ b/src/dynamics.f90 @@ -31,6 +31,7 @@ recursive subroutine FMS_Dynamics(Bundle, Timestep, GoalTime) real(kind=DefReal), intent(in) :: GoalTime, Timestep integer(kind=DefInt), parameter :: RecursionMax = 4 + integer(kind=DefInt) :: iTraj real(kind=DefReal) :: RedoGoalTime, RedoTimeStep logical, save :: FirstTime = .true. @@ -79,9 +80,32 @@ recursive subroutine FMS_Dynamics(Bundle, Timestep, GoalTime) !Stochastic collapse nuclear wavefunction !(RemoveDead must be called immediately after) call FMS_StochasticCollapse(Bundle) + write (fmiout, *) "-----------------------------------------------------------" + write (fmiout, *) "FMS_StochasticCollapse was successfully called, leaving now" + write (fmiout, *) "Now the time is: ", Bundle%CurrentTime + write (fmiout, *) "And these are our trajectories: (NumTraj, NCBFs)", Bundle%NumTraj, Bundle%NCBFs + write (fmiout, *) "And these are their centroids: " + do iTraj = 1, (((Bundle%NCBFs - 1) * Bundle%NCBFs) / 2) + write (fmiout, *) iTraj + write (fmiout, *) "is Cen to trajectories: ", Bundle%Centroids(iTraj)%CentID + write (fmiout, *) "Position", Bundle%Centroids(iTraj)%Particle(1)%get_pos() + end do + write (fmiout, *) "-----------------------------------------------------------" end if + write (fmiout, *) "Now calling RemoveDead... " call FMS_RemoveDead(Bundle) if (Bundle%NumTraj == 0) return + write (fmiout, *) "-----------------------------------------------------------" + write (fmiout, *) "Did deleting dead parts of Bundle mess up anything?" + write (fmiout, *) "Now the time is: ", Bundle%CurrentTime + write (fmiout, *) "And these are our trajectories: (NumTraj, NCBFs)", Bundle%NumTraj, Bundle%NCBFs + write (fmiout, *) "And these are their centroids: " + do iTraj = 1, (((Bundle%NCBFs - 1) * Bundle%NCBFs) / 2) + write (fmiout, *) iTraj + write (fmiout, *) "is Cen to trajectories: ", Bundle%Centroids(iTraj)%CentID + write (fmiout, *) "Position", Bundle%Centroids(iTraj)%Particle(1)%get_pos() + end do + write (fmiout, *) "-----------------------------------------------------------" end if ! Write output here only if user requested output at every step diff --git a/src/modules/BundleCalcsModule.f90 b/src/modules/BundleCalcsModule.f90 index 6b5dfaa..367a9d3 100644 --- a/src/modules/BundleCalcsModule.f90 +++ b/src/modules/BundleCalcsModule.f90 @@ -520,6 +520,13 @@ subroutine FMS_BuildHS(Bundle) real(kind=DefReal) :: time_tmp1, time_tmp2 call cpu_time(time_tmp1) + + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "We are building HS, check NumTraj, NCBFs" + !write (fmiout, *) "Number trajs:", Bundle%NumTraj + !write (fmiout, *) "Number CBFs:", Bundle%NCBFs + !write (fmiout, *) " ----------------------------------------------------" + Threshold = gldRegThresh ntraj = Bundle%NumTraj @@ -532,6 +539,15 @@ subroutine FMS_BuildHS(Bundle) H = (0., 0.) S_dot = (0., 0.) + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "Checking Centroid indices before calculating overlap: " + !write (fmiout, *) " ----------------------------------------------------" + !do i = 1, (((Bundle%NCBFs - 1) * Bundle%NCBFs) / 2) + ! write (fmiout, *) "CentID", Bundle%Centroids(i)%TrajID + ! write (fmiout, *) "Position", Bundle%Centroids(i)%Particle(1)%get_pos() + !end do + !write (fmiout, *) " ----------------------------------------------------" + ! Load H,S, and SDot matrices. do i = 1, ntraj do j = 1, i @@ -546,16 +562,18 @@ subroutine FMS_BuildHS(Bundle) ! TBF within the CBF. When the next Ms=2 TBF appears, we ! know that a new CBF is reached. if (glzCentroids .and. T_i%CBF /= T_j%CBF) then -! write(fmiOut,*) "i, j: ", i,j -! write(fmiOut,*) "IDi, IDj: ",T_i%StateID,T_j%StateID -! write(fmiOut,*) "Msi, Msj: ", T_i%Ms,T_j%Ms -! write(fmiOut,*) "CBFi, CBFj: ", T_i%CBF,T_j%CBF + write(fmiOut,*) "i, j: ", i,j + write(fmiOut,*) "IDi, IDj: ",T_i%StateID,T_j%StateID + write(fmiOut,*) "Msi, Msj: ", T_i%Ms,T_j%Ms + write(fmiOut,*) "CBFi, CBFj: ", T_i%CBF,T_j%CBF !! if( T_i%Ms.eq.2 .and. T_j%Ms.eq.2 ) then CBFi = T_i%CBF CBFj = T_j%CBF ICent = ((CBFi - 2) * (CBFi - 1)) / 2 + CBFj + write (fmiout, *) "This is what we are trying to use as CentID in BundleCalcs: ", ICent T_ij => Bundle%Centroids(ICent) + !write (fmiout, *) "Position", Bundle%Centroids(ICent)%Particle(1)%get_pos() !! endif end if ! GAIMS changed end @@ -574,31 +592,41 @@ subroutine FMS_BuildHS(Bundle) end if S(i, j) = overlap(T_i, T_j) + !write (fmiout, *) "I get here!!!" + !write (fmiout, *) "Are there Centroids?", glzCentroids ! Overlap will return elements of the individual overlap, which ! can be passed to kinetic. The overlaps themselves will be ! passed to Potential and to CGSDoTT instead of being recalculated. ! In case of Split-Operator: no need to calculate diagonal elements of H + if (glzFullyCoupled .or. IState /= JState) then ! If the overlap is beneath threshold, set matrix element to 0 if (abs(S(i, j)) >= FPZero) then + !write (fmiout, *) "----------------------------------------------------------------" if (i == j) then + !write (fmiout, *) "Pair of twice the same trajectory: i, j", i, j H(i, j) = overlap_KE(T_i, T_i, S(i, j)) + overlap_V(T_i, S_ij_precalc=S(i, j)) else if (glzCentroids) then + !write (fmiout, *) "Pair of different trajectories: i, j", i, j + !write (fmiout, *) "Calling overlap_V with existing centroids" + !write (fmiout, *) "This is the Centroid: ICent", ICent H(i, j) = overlap_KE(T_i, T_j, S(i, j)) + overlap_V(T_i, T_j, T_ij, S(i, j)) + !write (fmiout, *) "Got H(i, j)" else H(i, j) = overlap_KE(T_i, T_j, S(i, j)) + overlap_V(T_i, T_j, S_ij_precalc=S(i, j)) end if + !write (fmiout, *) "----------------------------------------------------------------" end if end if end if - ! fill in the offdiagonl + ! fill in the offdiagonal H(j, i) = conjg(H(i, j)) S(j, i) = conjg(S(i, j)) diff --git a/src/modules/OverlapModule.f90 b/src/modules/OverlapModule.f90 index 77cce87..51f6b20 100644 --- a/src/modules/OverlapModule.f90 +++ b/src/modules/OverlapModule.f90 @@ -1,4 +1,3 @@ -! Copyright Todd J. Martinez and Raphael D. Levine, 1994 module OverlapModule use GlobalModule use ParticleModule @@ -627,6 +626,9 @@ function overlap_V_trajectory(T_i, T_j, T_c, S_ij_precalc) result(V) !bfec if (glzCentroids .and. (CBFi /= CBFj)) then ic = T_c%CentID(1); jc = T_c%CentID(2) + write (fmiout, *) " ------------------------------------------------------------ " + write (fmiout, *) "In the OverlapModule: getting ic and jc", ic, jc + write (fmiout, *) " ------------------------------------------------------------ " end if is = T_i%StateID; js = T_j%StateID @@ -672,7 +674,8 @@ function overlap_V_trajectory(T_i, T_j, T_c, S_ij_precalc) result(V) ! check that the correct centroid was passed ! cent_match = (it==ic .and. jt==jc) .or. (jt==ic .and. it==jc) cent_match = (CBFi == ic .and. CBFj == jc) .or. (CBFj == ic .and. CBFi == jc) - write (fmiout, *) "(CBFij == ijc .and. CBFji == jic)", "CBFi, ic, CBFj, jc", CBFi, ic, CBFj, jc + !write (fmiout, *) "(CBFij == ijc .and. CBFji == jic)", "CBFi, ic, CBFj, jc", CBFi, ic, CBFj, jc + !write (fmiout, *) "Which trajectories is T_ij (=T_c) centroid to?", T_c%CentID(1), T_c%CentID(2) flush (fmiout) ! cent_match = (i==ic .and. j==jc) .or. (j==ic .and. i==jc) if (.not. cent_match) then diff --git a/src/modules/SamplingModule.f90 b/src/modules/SamplingModule.f90 index 79c9be3..5946f68 100644 --- a/src/modules/SamplingModule.f90 +++ b/src/modules/SamplingModule.f90 @@ -1341,6 +1341,8 @@ subroutine FMS_InitialSwarm(B1) real(defReal) :: B1_norm complex(DefComp), dimension(B1%NumTraj) :: S_if ! overlap between initial and final + write (fmiout, *) "Initializing a swarm of trajectories" + ntraj = B1%NumTraj natom = B1%NumParticles nstate = B1%NumStates @@ -1356,18 +1358,42 @@ subroutine FMS_InitialSwarm(B1) P_norm = sqrt(sum(P**2)) ! 2. set up the Bundle trajectories - do n = 1, ntraj - call random_number(dX) - dX = sigma * (2.0d0 * dX - 1.d0) + ! -- Added just for usage with SOC model and sepSS testing: --------------------------------------- + if (gliModel == FMSZERO) then + write (fmiout, *) "Initialising for ToyModel in 1D" + dX = [0.5d0, 0.d0, 0.d0] + dP = [0.d0, 0.d0, 0.d0] + write (fmiout, *) "dX, dP: ", dX, dP + + if (ntraj > 2) then + do n = 2, ntraj + B1%Trajectory(n)%TrajID = n + call B1%Trajectory(n)%set_pos(X + dX) + call B1%Trajectory(n)%set_mom(P + dP) ! useless here because dP is zero, but for completeness... + dX = [10.0d0, 0.d0, 0.d0] + end do + else + call B1%Trajectory(ntraj)%set_pos(X + dX) + call B1%Trajectory(ntraj)%set_mom(P + dP) ! useless here because dP is zero, but for completeness... + end if - call random_number(dP) - dP = sigma / P_norm * (2.0d0 * dP - 1.d0) + ! ------------------------------------------------------------------------------------------------- - B1%Trajectory(n)%TrajID = n - call B1%Trajectory(n)%set_pos(X + dX) - call B1%Trajectory(n)%set_mom(P + dP) - end do + else + do n = 1, ntraj + + call random_number(dX) + dX = sigma * (2.0d0 * dX - 1.d0) + + call random_number(dP) + dP = sigma / P_norm * (2.0d0 * dP - 1.d0) + + B1%Trajectory(n)%TrajID = n + call B1%Trajectory(n)%set_pos(X + dX) + call B1%Trajectory(n)%set_mom(P + dP) + end do + end if ! 3. update centroids !bfec @@ -1380,36 +1406,44 @@ subroutine FMS_InitialSwarm(B1) end do end if - ! normalize the Bundle - B1_norm = FMS_Norm(B1) - ! TODO: This should be a bundle method - do n = 1, ntraj - B1%Trajectory(n)%Amplitude = B1%Trajectory(n)%Amplitude / sqrt(B1_norm) + do n = 1, B1%NumTraj + write (fmiout, *) "Before normalising Bundle: Amplitude of traj ", n, " :", B1%Trajectory(n)%Amplitude end do - ! 3. work out the overlaps between the Bundle trajectorirs - call FMS_BuildHS(B1) + !! normalize the Bundle + !B1_norm = FMS_Norm(B1) + !! TODO: This should be a bundle method + !do n = 1, ntraj + ! B1%Trajectory(n)%Amplitude = B1%Trajectory(n)%Amplitude / sqrt(B1_norm) + !end do - ! 4. overlap the intial Gaussian from Geometry.dat onto the Bundle - do n = 1, ntraj - S_if(n) = overlap(T_init, B1%Trajectory(n)) - end do + !! 3. work out the overlaps between the Bundle trajectories + !call FMS_BuildHS(B1) + + !! 4. overlap the intial Gaussian from Geometry.dat onto the Bundle + !do n = 1, ntraj + ! S_if(n) = overlap(T_init, B1%Trajectory(n)) + !end do - ! careful here, matmul( A, B ) = A^* . B - S_if = matmul(conjg(FMS_bSInvMat(B1)), S_if) + !! careful here, matmul( A, B ) = A^* . B + !S_if = matmul(conjg(FMS_bSInvMat(B1)), S_if) - ! 5. set the amplitudes - call FMS_Set_Amplitude(B1, S_if) + !! 5. set the amplitudes + !call FMS_Set_Amplitude(B1, S_if) call T_init%destroy() B1_norm = FMS_Norm(B1) - write (fmiOut, *) "Overlap between Initial Trajectory and the Bundle" - write (fmiOut, *) sqrt(B1_norm) + !write (fmiOut, *) "Overlap between Initial Trajectory and the Bundle" + !write (fmiOut, *) sqrt(B1_norm) ! normalize the Bundle again do n = 1, ntraj B1%Trajectory(n)%Amplitude = B1%Trajectory(n)%Amplitude / sqrt(B1_norm) end do + + do n = 1, B1%NumTraj + write (fmiout, *) "After normalising Bundle: Amplitude of traj ", n, " :", B1%Trajectory(n)%Amplitude + end do end subroutine FMS_InitialSwarm end module SamplingModule diff --git a/src/modules/SelectionModule.f90 b/src/modules/SelectionModule.f90 index b155f7a..031d370 100644 --- a/src/modules/SelectionModule.f90 +++ b/src/modules/SelectionModule.f90 @@ -59,6 +59,25 @@ subroutine FMS_StochasticCollapse(B1) type(T_TrajectoryBundle), intent(inout) :: B1 ! An array of temporary Bundles for single state stochastic selection type(T_TrajectoryBundle), allocatable :: BundleSS(:) +! sepSS: array of temporary Bundles for singlet and triplet separated stochastic selection +! and the dimension of sep bundles array + type(T_TrajectoryBundle), allocatable :: BundlesMult(:) + integer(kind=DefInt) :: BundlesMultDim +! --- Added 4 sep SS +! sepSS: needed to determine total Coupled matrix before stochastic selection is performed +! (same as in FMS_BuildCoupled and FMS_GroupIntoBlocks which are otherwise +! only called from perform_stochastic_selection) + integer(kind=DefInt), dimension(B1%NumTraj, B1%NumTraj) :: Coupled + integer(kind=DefInt), dimension(B1%NumTraj + 1, B1%NumTraj + 1) :: blocktrajid + integer(kind=DefInt), dimension(B1%NumTraj + 1, B1%NumTraj + 1) :: MergedBlockTrajID, MergedSepBlockTrajID + integer(kind=DefInt), dimension(B1%NumTraj + 1) :: ntrajblock + integer(kind=DefInt) :: nblockMerged, DeadID + integer(kind=DefInt), dimension(50) :: SaveForceKill +! sepSS: Shall we perform multiplet separated stochastic selection + integer(kind=DefInt), allocatable :: StoSelMode(:), NumTrajBlock(:) + logical :: isSing, isTrip, isSingTrip +! --- Added 4 sep SS end --- + integer(kind=DefInt) :: iTraj, i, j, nblock ! AIMSWISS: Shall we perform stochastic selection logical :: performSelection ! AIMSWISS: Current selection time @@ -106,6 +125,254 @@ subroutine FMS_StochasticCollapse(B1) end do deallocate (BundleSS) +! Multiplicity-Specific Stochastic Selection +! 1. Decide whether we have Sing only, Trip only, or Sing and Trip +! 2. Create temporary, separate Bundles that only hold the trajectories within a multiplicity +! 3. Perform selections within those bundles +! 4. Reunite to reform the original bundle, get rid of temporary block Bundles + + else if (NTrip /= 0) then + + ! (1) getting Coupled matrix to identify separated S only, or T only blocks + + ! - get the total converged coupled matrix, and the number of blocks + write (fmiout, *) "-----------------------------------------------------------" + write (fmiout, *) "Building coupled matrix for total Bundle (B1, unseparated)" + Coupled = 0 + !write (fmiout, *) "This is the total, converged coupled matrix of size ", size(Coupled) + call FMS_BuildCoupled(B1, Coupled, selectionTime) + write (fmiout, *) "This is the total, converged coupled matrix of size ", size(Coupled) + do iTraj = 1, B1%NumTraj + write (fmiout, *) Coupled(iTraj,:) + end do + + write (fmiout, *) "Get number of blocks for total Bundle (B1, unseparated)" + blocktrajid = 0 + ntrajblock = 0 + nblock = 0 + call FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, ntrajblock, nblock) + write (fmiout, *) "Number of Blocks: ", nblock + write (fmiout, *) "-----------------------------------------------------------" + + ! - merge the total coupled matrix according to + ! - Singlet only --> perform StoSel + ! - Triplet only --> perform StoSel + ! - Singlet and Triplet in the same block --> no StoSel on these mixed blocks + isSing = .false. + isTrip = .false. + isSingTrip = .false. + + write (fmiout, *) "-----------------------------------------------------------" + write (fmiout, *) "Deciding StoSel mode based on total Bundle" + + allocate (StoSelMode(nblock)) + allocate (NumTrajBlock(nblock)) + + ! Get isSIng, isTrip, isSingTrip for each block in the original, total Coupled matrix + NumTrajBlock = 0 + do i = 1, nblock + isSing = .false. + isTrip = .false. + isSingTrip = .false. + do iTraj = 1, B1%NumTraj + if (any(blocktrajid(:,i) == iTraj)) then + NumTrajBlock(i) = NumTrajBlock(i) + 1 + if (B1%Trajectory(iTraj)%triplet) then + isTrip = .true. + else + isSing = .true. + end if + + end if + end do + if (isTrip .eqv. .true. .and. isSing .eqv. .true.) then + isSingTrip = .true. + StoSelMode(i) = 0 ! no selection because Sing and Trip are coupled + else if (isSing .eqv. .true.) then + StoSelMode(i) = 1 ! selection for Sing + else if (isTrip .eqv. .true.) then + StoSelMode(i) = 3 ! selection for Trip + else + call FMS_DieError("No valid Selection mode could be determined") + end if + + write (fmiout, *) "For block ", i, "StoSel mode is: " + write (fmiout, *) "isSing is ", isSing, ", isTrip is ", isTrip, ", and... isSingTrip is", isSingTrip + write (fmiout, *) "This is for the trajs ", blocktrajid(:,i) + write (fmiout, *) "This is the current StoSel mode: ", StoSelMode(i) + end do + + ! Merge blocks depending on isSIng, isTrip + ! (No action for isSingTrip because not to be touched in StoSel + do i = 1, nblock + + if (StoSelMode(i) == 0) cycle + + do j = i+1, nblock + + if (StoSelMode(i) == StoSelMode(j)) then + ! Copy over block j TrajIDs to i column of blocktrajid + blocktrajid(NumTrajBlock(i)+1:NumTrajBlock(i)+NumTrajBlock(j), i) = blocktrajid(1:NumTrajBlock(j),j) + ! Set j column to zero + blocktrajid(:,j) = 0 + ! Set StoSelMode to zero for block copied over (the j column) + StoSelMode(j) = 0 + + end if + + end do + + end do + MergedBlockTrajID = blocktrajid + nblockMerged = nblock + ! Generate 2nd version of MergedBlockTrajID w/ indices ref to what will become the sep bundles + ! This is the MergedSepBlockTrajID matrx (sorry complicated.. :( ) + ! TODO: this doesn't seem robust at all, go back later and improve + MergedSepBlockTrajID = MergedBlockTrajID + do i = 1, B1%NumTraj + 1 + do j = 2, B1%NumTraj + 1 + if (MergedBlockTrajID(i,j) == 0) cycle + MergedSepBlockTrajID(i,j) = MergedBlockTrajID(i,j) - MergedBlockTrajID(1,j) + 1 + end do + end do + + write (fmiout, *) "This is the S, T, S/T StoSel mode merged blocktrajid: " + do j = 1, B1%NumTraj + 1 + write (fmiout, *) blocktrajid(j,:) + end do + write (fmiout, *) "And these are our StoSel (updated) modes: ", StoSelMode + + write (fmiout, *) "This is the merged sepblocktrajid with what will be sep bundle IDs: " + do j = 1, B1%NumTraj + 1 + write (fmiout, *) MergedSepBlockTrajID(j,:) + end do + write (fmiout, *) "And these are our StoSel (updated) modes: ", StoSelMode + + ! (2) Create and fill temp separated bundles + + ! Determine how many block bundles we need to create, based on merged blocktrajid: + ! and based on whether we need to carry out SS in the merged block or not + BundlesMultDim = 0 + do i = 1, nblock + if (StoSelMode(i) /= 0) then + BundlesMultDim = BundlesMultDim + 1 + end if + end do + nblock = BundlesMultDim + write (fmiout, *) "BundlesMultDim based on merged blocktrajid: ", nblock + write (fmiout, *) "nblock based on merged blocktrajid: ", nblock + + write (fmiout, *) "-----------------------------------------------------------" + + allocate (BundlesMult(BundlesMultDim)) + + ! Create BundlesMultDim many temp bundles and fill them up + ! with corresponding trajs from B1 + call getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) + + ! (3) Perform SS on the separated bundles + do i = 1, nblock + + if (StoSelMode(i) == 1 .or. StoSelMode(i) == 3) then + + write (fmiout, *) "StoSelMode for block ", i, "is ", StoSelMode(i) + ! Do the stochastic selection + if (BundlesMult(i)%NCBFs > 1 ) then + write (fmiout, *) "Performing stochastic selection for block Bundle", i, & + "with", BundlesMult(i)%NCBFs, "CBFs" + + call perform_stochastic_selection(BundlesMult(i), selectionTime) ! TODO: make sure the right Trajs + ! get killed (currently: gliForceKill + ! uses TrajID based on sep bundle 2 + ! kill in B1 + + if (i>1) then + write (fmiout, *) "Adjusting gliForceKill using MergedSepBlockTrajID (this is the full matrix): " + do j = 1, B1%NumTraj + 1 + write (fmiout, *) MergedSepBlockTrajID(j,:) + end do + + write (fmiout, *) "and MergedBlockTrajID (this is the full matrix): " + do j = 1, B1%NumTraj + 1 + write (fmiout, *) MergedBlockTrajID(j,:) + end do + + write (fmiout, *) "gliForcKill for block Bundle", i + write (fmiout, *) gliForceKill + write (fmiout, *) "adjusting to B1 indiced..." + + SaveForceKill = gliForceKill + + do iTraj = 2, B1%NumTraj + 1 + do j = 1, B1%NumTraj + 1 + if (MergedSepBlockTrajID(iTraj,j) == 0) cycle + do DeadID = 1, size(SaveForceKill) + write (fmiout, *) "many numbers (?)", SaveForceKill(DeadID) + if (MergedSepBlockTrajID(iTraj,j) == SaveForceKill(DeadID)) then + write (fmiout, *) "Match for gliForceKill: ", SaveForceKill(DeadID), "(this should be 4)" + write (fmiout, *) "with MergedSepblock: ", MergedSepBlockTrajID(iTraj,j) + write (fmiout, *) "i and j are: ", iTraj, j, "(should be 2, 4)" + write (fmiout, *) "Corresponding MergedBlock: ", MergedBlockTrajID(iTraj,j), "(this should be 6)" + gliForceKill(DeadID) = MergedBlockTrajID(iTraj,j) + end if + end do + end do + end do + + write (fmiout, *) "Adjusted gliForcKill for block Bundle", i, "and obtained" + write (fmiout, *) gliForceKill + end if + + else + write (fmiout, *) "No stochastic selection because block Bundle only has one CBF" + end if + + else + if (any(blocktrajid(:,i) /= 0)) then + write (fmiout, *) "No stochastic selection because block Bundle", i, "contains a mix of S and T" + else + write (fmiout, *) "No stochastic selection because block Bundle", i, "was merged into other block" + end if + end if + + end do + + !write (fmiout, *) "Check that Centroids of B1 remain unaffected" + !write (fmiout, *) "Part 1: B1 Centroids before" + !do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) + ! write (fmiout, *) "Centroid number", iTraj + ! write (fmiout, *) "Is centroid to trajectories", B1%Centroids(iTraj)%CentID + ! write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() + !end do + + ! (4) Set appropriate amplitudes to zero in the original bundle + ! If selection happened, copy information which trajs died + if (any(gliForceKill(:) /= 0)) then + call copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode) + write (fmiout, *) "States were copied back without error" + end if + + !write (fmiout, *) "Check that Centroids of B1 remain unaffected" + !write (fmiout, *) "Part 2: B1 Centroids after" + !do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) + ! write (fmiout, *) "Centroid number", iTraj + ! write (fmiout, *) "Is centroid to trajectories", B1%Centroids(iTraj)%CentID + ! write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() + !end do + + ! Destroy the temporary block bundles + do i = 1, nblock + call BundlesMult(i)%destroy() + end do + + !call FMS_DieError("Now, please go back and reconsider your life choices!") + + !else + + ! call perform_stochastic_selection(B1, selectionTime) + + !end if + else call perform_stochastic_selection(B1, selectionTime) @@ -127,7 +394,7 @@ subroutine perform_stochastic_selection(B1, selectionTime) ! Coupling matrix for TBF basis integer(kind=DefInt), dimension(B1%NumTraj, B1%NumTraj) :: Coupled ! Number of blocks - integer(kind=DefInt) :: nblock + integer(kind=DefInt) :: nblock, iTraj ! Matrix containing IDs of TBFs belonging to a block. integer(kind=DefInt), dimension(B1%NumTraj + 1, B1%NumTraj + 1) :: blocktrajid ! Array containing number of TBFs belonging to a block. @@ -141,6 +408,13 @@ subroutine perform_stochastic_selection(B1, selectionTime) ! AIMSWISS: Current selection time real(kind=DefReal), intent(in) :: selectionTime + write (fmiout, *) "Number of traj used for dimension of Coupled: ", B1%NumTraj + write (fmiout, *) "meaning Coupled dimension is: ", "( ", B1%NumTraj, ", ", B1%NumTraj, " )" + + write (fmiout, *) "CentIDs in perform_stochastic_selection (iTraj, iCBF:" + do iTraj = 1, B1%NumTraj + write (fmiout, *) B1%Trajectory(iTraj)%TrajID, B1%Trajectory(iTraj)%CBF + end do ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ! First, decompose the hamiltonian into a block diagonal representation ! work out the coupling matrix @@ -153,6 +427,7 @@ subroutine perform_stochastic_selection(B1, selectionTime) ntrajblock = 0 nblock = 0 call FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, ntrajblock, nblock) + write (fmiout, *) "Number of Blocks: ", nblock ! There should be at least one block of trajectories if (nblock == 0) then @@ -185,6 +460,7 @@ subroutine perform_stochastic_selection(B1, selectionTime) ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ! Remove all TBFs that do not belong selected block call FMS_KillOtherBlocks(B1, nblock, ntrajblock, iblockslct, blocktrajid) + write (fmiout, *) "Leaving perform_stochastic_selection" end subroutine perform_stochastic_selection @@ -341,6 +617,10 @@ subroutine FMS_BuildCoupled(B1, Coupled, selectionTime) real(kind=DefReal), intent(in) :: selectionTime integer(kind=DefInt) :: ntraj, i, j + write (fmiout, *) "---------------------------------------------------------------------------" + write (fmiout, *) " Inside FMS_BuildCoupled " + write (fmiout, *) "The number of traj we are using is: ", B1%NumTraj + ntraj = B1%NumTraj do i = 1, ntraj @@ -350,6 +630,8 @@ subroutine FMS_BuildCoupled(B1, Coupled, selectionTime) if (B1%Trajectory(i)%cbf == B1%Trajectory(j)%cbf) then Coupled(i, j) = 1 Coupled(j, i) = 1 + write (fmiout, *) "Coupled(i,j) set to 1 because trajs have the same CBF" + write (fmiout, *) "CBFi, CBFj", B1%Trajectory(i)%cbf, B1%Trajectory(j)%cbf cycle end if end do @@ -365,9 +647,17 @@ subroutine FMS_BuildCoupled(B1, Coupled, selectionTime) else + do i = 1, ntraj + write (fmiout, *) "Building coupled for a triplet?", B1%Trajectory(i)%triplet + end do + call FMS_BuildCoupled_ESS(B1, Coupled) + write (fmiOut, *) "Coupled matrix is:", Coupled, ", Size is:", size(Coupled) end if + + write (fmiout, *) "---------------------------------------------------------------------------" + end subroutine FMS_BuildCoupled ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -380,8 +670,18 @@ subroutine FMS_BuildCoupled_ESS(B1, Coupled) ntraj = B1%NumTraj + write (fmiout, *) "Calculating coupling for trajs in bundle" + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "Are traj indices still fine here?" + write (fmiout, *) "Number trajs:", B1%NumTraj + write (fmiout, *) "Number CBFs:", B1%NCBFs + write (fmiout, *) " ----------------------------------------------------" + do i = 2, ntraj do j = 1, i + write (fmiout, *) "ID, CBF, is triplet?" + write (fmiout, *) "i", B1%Trajectory(i)%TrajID, B1%Trajectory(i)%CBF, B1%Trajectory(i)%triplet + write (fmiout, *) "j", B1%Trajectory(j)%TrajID, B1%Trajectory(j)%CBF, B1%Trajectory(j)%triplet if (B1%Trajectory(i)%cbf == B1%Trajectory(j)%cbf) cycle if (abs(FMS_bH(B1, i, j)) > gldStochaThresh) then Coupled(i, j) = 1 @@ -391,6 +691,7 @@ subroutine FMS_BuildCoupled_ESS(B1, Coupled) end do ! Uncover all indirect connections between TBFs +!!!!!!!!!!!!!!! Look at thattt: call FMS_ConvergeCoupled(ntraj, Coupled) end subroutine FMS_BuildCoupled_ESS @@ -514,6 +815,9 @@ subroutine FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, & integer(kind=DefInt) :: ntraj, i, j, l, itraj + write (fmiout, *) "=========================================" + write (fmiout, *) "What happens in GroupIntoBlocks?" + ntraj = B1%NumTraj ztrajdone(:) = .false. @@ -527,22 +831,31 @@ subroutine FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, & if (ztrajdone(j)) cycle blocktrajid(1, i) = j ztrajdone(j) = .true. + write (fmiout, *) "The hitherto unsorted traj ", j, "was added to the block as the first traj" exit end do itraj = blocktrajid(1, i) !first trajectory in block + write (fmiout, *) "This is blocktrajid inside the i loop, i = ", i + do j = 1, B1%NumTraj + 1 + write (fmiout, *) blocktrajid(j,:) + end do + ! If no trajectories added, then exit loop if (itraj == 0) then nblock = i - 1 + write (fmiout, *) "nblock is now ", nblock exit end if ! Add other coupled trajectories to this current block l = 1 do j = 1, ntraj + write (fmiout, *) "Now sorting traj j = ", j if (ztrajdone(j)) cycle !if this trajectory has already been sorted skip if (Coupled(j, itraj) == 1) then + write (fmiout, *) "traj ", j, "is in the same block as traj ", itraj l = l + 1 blocktrajid(l, i) = j ztrajdone(j) = .true. @@ -552,6 +865,16 @@ subroutine FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, & end do + write (fmiout, *) "This is the finished blocktrajid: " + do j = 1, B1%NumTraj + 1 + write (fmiout, *) blocktrajid(j,:) + end do + + write (fmiout, *) "These trajectories are in the first block: ", blocktrajid(:,1) + write (fmiout, *) "and these are in the second block: ", blocktrajid(:,2) + + write (fmiout, *) "=========================================" + end subroutine FMS_GroupIntoBlocks ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -692,19 +1015,30 @@ subroutine FMS_KillOtherBlocks(B1, nblock, ntrajblock, iblockslct, blocktrajid) integer(kind=DefInt) :: i, j, l, jtraj ! Mark all other trajectories for death l = 0 + write (fmiout, *) "-----------------------------------------------" + write (fmiout, *) "Marking trajectories for death" do i = 1, nblock + write (fmiout, *) "Current block: ", i + write (fmiout, *) "was it selected? ", i == iblockslct if (i == iblockslct) cycle do j = 1, ntrajblock(i) + write (fmiout, *) "We have a j! ", j l = l + 1 jtraj = blocktrajid(j, i) gliForceKill(l) = B1%Trajectory(jtraj)%TrajID + write (fmiout, *) "Trajectory will be killed (TrajID, CBF): ", B1%Trajectory(jtraj)%TrajID, B1%Trajectory(jtraj)%CBF B1%Trajectory(jtraj)%Amplitude = dcmplx(0.0d0, 0.0d0) end do end do + write (fmiout, *) "Array of TrajIDs for dead/ non dead Trajs: ", gliForceKill, "size is: ", size(gliForceKill) + write (fmiout, *) "-----------------------------------------------" ! Remove any dead trajectory amplitudes do jtraj = 1, B1%NumDeadTraj B1%DeadTraj(jtraj)%Amplitude = dcmplx(0.0d0, 0.0d0) + write (fmiout, *) "Killing this trajectory: ", jtraj + write (fmiout, *) "Out of this many total dead trajs: ", B1%NumDeadTraj + write (fmiout, *) "Was it a triplet?", B1%DeadTraj(jtraj)%triplet end do end subroutine FMS_KillOtherBlocks @@ -764,6 +1098,400 @@ subroutine FMS_WriteSelectionLog(parent, child, currentS, predictedS, l_select) end subroutine FMS_WriteSelectionLog + subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) + + ! Number of separated bundles we need to create + ! is determined based on merged blocktrajid + ! so this number is also how many Selections we will perform + integer(kind=DefInt), intent(in) :: BundlesMultDim + ! the original bundle + type(T_TrajectoryBundle), intent(in) :: B1 + ! 2D array holding the TrajIDs sorted according to StoSel blocks they belong to + integer(kind=DefInt), intent(in), dimension(B1%NumTraj + 1, B1%NumTraj + 1) :: blocktrajid + ! array holding the separated bundles + type(T_TrajectoryBundle), intent(inout) :: BundlesMult(BundlesMultDim) + + ! Getting info on Traj, CBF, Cen indices for each StoSel block: + ! - Traj information + integer(kind=DefInt), dimension(BundlesMultDim) :: NumTrajBlock + integer(kind=DefInt) :: iTraj, addTraj + ! - CBF information + integer(kind=DefInt), dimension(B1%NCBFs + 1, BundlesMultDim) :: BlockCBFID + integer(kind=DefInt), dimension(BundlesMultDim) :: NumCBFBlock + integer(kind=DefInt) :: iblocks, jblocks, BlockCBFCount + integer(kind=DefInt) :: CBFcount, CBFcurr, CBF_i, CBF_j + ! - Cen information + integer(kind=DefInt) :: Cen_i, Cen_j, CenCount, Cen_TrajID, CenID + logical :: is_addCen + ! Number of particles + integer(kind=DefInt) :: npart + + npart = B1%Trajectory(1)%NumParticles + + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "Total number of trajs and CBFs in unseparated Bundle: " + write (fmiout, *) B1%NumTraj, B1%NCBFs + write (fmiout, *) " ----------------------------------------------------" + + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "Number of trajs and CBFs that will be in unseparated Bundle: " + !write (fmiout, *) "Number S, Number T trajs:", NumSingTraj, NumTripTraj + !write (fmiout, *) "Number CBFs S, Number CBFs T trajs:", nSingCBF, nTripCBF + !write (fmiout, *) " ----------------------------------------------------" + + ! Determine number of trajs and CBFs in each StoSel block + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "What we currently have in the B1 bundle: " + do iTraj = 1, B1%NumTraj + + write (fmiout, *) "Info for trajectory with ID ", B1%Trajectory(iTraj)%TrajID + write (fmiout, *) "Triplet?", B1%Trajectory(iTraj)%triplet + write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms + write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF + + end do + + NumTrajBlock = 0 + NumCBFBlock = 0 + BlockCBFID = 0 + CBFcurr = 0 + CBFcount = 0 + + ! Count how many trajs and CBFs we have in the current block + ! - Construct NumTrajBlock (will be needed 4 allocating temp Bundle and + ! knowing how many trajs 2 add to temp Bundles) + ! - Construct BlockCBFID 2 keep track of CBFs corresponding 2 trajs + ! and NumCBFBlock (will be needed 4 allocating temp Bundle and + ! resetting CBF IDs in these temp Bundles) + ! - CBFcurr, CBFcount used to make sure that multiplets are copied + ! over correctly, and not separated + do iblocks = 1, bundlesmultdim + + CBFcount = 0 + do iTraj =1, B1%NumTraj + + ! count traj if its ID is in the current block + if (any(blocktrajid(:, iblocks) == iTraj)) then + NumTrajBlock(iblocks) = NumTrajBlock(iblocks) + 1 + if (B1%Trajectory(iTraj)%CBF /= CBFcurr) then + CBFcurr = B1%Trajectory(iTraj)%CBF + CBFcount = CBFcount + 1 + BlockCBFID(CBFcount, iblocks) = CBFcurr + NumCBFBlock(iblocks) = NumCBFBlock(iblocks) + 1 + end if + end if + + end do + + write (fmiout, *) "Number of Trajs and CBFs in block ", iblocks + write (fmiout, *) NumTrajBlock(iblocks), NumCBFBlock(iblocks) + + end do + + write (fmiout, *) "This is BlockCBFID :" + do iblocks = 1, BundlesMultDim + write (fmiout, *) BlockCBFID(:, iblocks) + end do + + !Create temporary bundles, separated according 2 selection mode: + ! S only, or T only, (as for S/T combined, no SS should happen + ! no temp Bundles are created 4 this case) + do iblocks = 1, BundlesMultDim + + call BundlesMult(iblocks)%create(numtraj=NumTrajBlock(iblocks), & + numdeadtraj=0, & !!!!! TODO: Is it fine setting this to 0??? + numstates=B1%NumStates, & + numparticles=npart, & + ncbfs=NumCBFBlock(iblocks)) + + end do + + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "These are the empty, separated Bundles we created: " + !do iblocks = 1, BundlesMultDim + ! write (fmiout, *) "Block number ", iblocks + ! write (fmiout, *) "Number trajs, Number CBFs:", BundlesMult(iblocks)%NumTraj, BundlesMult(iblocks)%NCBFs + !end do + !write (fmiout, *) " ----------------------------------------------------" + + ! Copy over to the temp, separated bundles: + ! - matching trajs according 2 BlockTrajID column + ! - matching Centroids (make sure CentID matches) + write (fmiout, *) "The blocktrajid we are using for copying trajs to sep Bundles: " + do iTraj = 1, B1%NumTraj + 1 + write (fmiout, *) blocktrajid(iTraj,:) + end do + + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "Original centroids in B1" + write (fmiout, *) "Does B1 have smth like positions for Centroids?", B1%Centroids(1)%Particle(1)%get_pos() + do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) + write (fmiout, *) "Centroid number ", iTraj + write (fmiout, *) "is centroid to trajectories ", B1%Centroids(iTraj)%CentID + write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() + end do + + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "Filling up the empty, separated bundles with Trajs from B1:" + + do iblocks = 1, BundlesMultDim + + write (fmiout, *) "=====================================================" + write (fmiout, *) "Now adding trajs to sep Bundle ", iblocks + write (fmiout, *) "Should start with traj of index ", blocktrajid(1,iblocks) + write (fmiout, *) "=====================================================" + + CenCount = 1 + CBFcount = 1 + + do iTraj = 1, NumTrajBlock(iblocks) + + addTraj = blocktrajid(iTraj,iblocks) + ! iTraj --> ID traj will have in temp bundle + ! addTraj --> ID of corresponding traj in B1 + + write (fmiout, *) "iTraj ", iTraj + write (fmiout, *) "represents traj", addTraj + write (fmiout, *) "The BlockCBFID for current block: ", BlockCBFID(:, iblocks) + write (fmiout, *) "has actual CBF", B1%Trajectory(addTraj)%CBF, "the BlockCBFID CBF is:", BlockCBFID(CBFcount,iblocks) + + !if (iblocks>1) then + ! write (fmiout, *) "in sep Bundle, needs to have CBF", B1%Trajectory(addTraj)%CBF - NumCBFBlock(iblocks-1) + !else + ! write (fmiout, *) "in sep Bundle, needs to have CBF", B1%Trajectory(addTraj)%CBF + !end if + + !write (fmiout, *) "We are counting CBF (CBFcount) ", CBFcount + !write (fmiout, *) B1%Trajectory(iTraj)%CBF, BlockCBFID(CBFcount,iblocks) + + if (addTraj == 0) cycle ! dont try to add traj if there is no corresponding traj in BlockTrajID + ! TODO: do we actually need this? + + ! increase CBFcount if CBFcurr has changed from previous iteration + ! (BlockCBFID(CBFcount,iblocks) is the CBF of previous the iteration) + if (B1%Trajectory(addTraj)%CBF /= BlockCBFID(CBFcount,iblocks)) then + !CBFcurr = B1%Trajectory(iTraj)%CBF + write (fmiout, *) "Adding one to CBFcount" + CBFcount = CBFcount + 1 + end if + + ! create the traj in the temp separated bundle and set it 2 corresponding traj in B1 + ! set TrajID to iTraj to make sure indices in temp bundle start at 1 + call BundlesMult(iblocks)%Trajectory(iTraj)%create(npart, B1%NumStates) + BundlesMult(iblocks)%Trajectory(iTraj) = B1%Trajectory(addTraj) + BundlesMult(iblocks)%Trajectory(iTraj)%TrajID = iTraj + + ! check if CBF IDs need to be reset to 1 (the case if we are not in the first separated block) + ! to make sure CBF IDs also start at 1 + if (iblocks>1) then + BlockCBFCount = 0 + do jblocks = 1, iblocks-1 + BlockCBFCount = BlockCBFCount + NumCBFBlock(jblocks) + end do + write (fmiout, *) "BlockCBFcount (all CBFs added to previous blocks: ", BlockCBFCount + BundlesMult(iblocks)%Trajectory(iTraj)%CBF = B1%Trajectory(addTraj)%CBF - BlockCBFCount + write (fmiout, *) "CBF receives the number (must start at 1)", B1%Trajectory(addTraj)%CBF - BlockCBFCount + else + BundlesMult(iblocks)%Trajectory(iTraj)%CBF = B1%Trajectory(addTraj)%CBF + end if + + ! Next, copy over matching centroids from B1 to temp bundle + ! --------------------------------------------------------------------------------------- + ! Careful here, figuring this out has been a mess but we (we=Vera) pray that it works now + ! - do we really cover all possible CBF ID pairs, i.e. catch all centroids? + ! --------------------------------------------------------------------------------------- + write (fmiout, *) "___________________________" + + ! - start from the 1st CBF ID (i.e. 1st row, current column of BlockCBFID) + ! and go over all CBF IDs as CBFcount increases + CBF_i = BlockCBFID(CBFcount, iblocks) + write (fmiout, *) "CBF_i is, CBF_j will be", CBF_i, CBF_i + 1 + write (fmiout, *) "loop limit will be:", B1%NCBFs + + ! - loop over possible CBF ID pairs and check 4 matching CentID pair in B1 + is_addCen = .false. + do CBF_j = CBF_i + 1, B1%NCBFs ! what if CBF IDs are not in ascending order, does this still work??? :S + + if (CenCount > (((BundlesMult(iblocks)%NCBFs - 1) * BundlesMult(iblocks)%NCBFs) / 2) ) cycle ! dont add more Cens + ! if max of cens in + ! temp bundle is reached + + write (fmiout, *) "CBF_i and CBF_j", CBF_i, CBF_j + if (any(CBF_i == BlockCBFID(:,iblocks)) .and. any(CBF_j == BlockCBFID(:,iblocks)) .and. CBF_i /= CBF_j) then + + write (fmiout, *) "dream come true", CBF_j, CBF_i + write (fmiout, *) "CenID is", CenID + + ! For this CentID pair, find out corresponding ID of Centroid in B1 (CenID) + do Cen_TrajID = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) + + write(fmiout, *) "Is Cen ", Cen_TrajID, "the correct one?" + write (fmiout, *) B1%Centroids(Cen_TrajID)%CentID + Cen_i = B1%Centroids(Cen_TrajID)%CentID(1); Cen_j = B1%Centroids(Cen_TrajID)%CentID(2) + + if (CBF_i == Cen_i .and. CBF_j == Cen_j .or. CBF_j == Cen_i .and. CBF_i == Cen_j) then + CenID = Cen_TrajID + write (fmiout, *) "The ID of the centroid we need is ", CenID + is_addCen = .true. + !cycle + end if + + end do + + if (is_addCen .eqv. .true.) then + ! Use the corresponding ID of Centroid in B1 (CenID) to copy over matching info + BundlesMult(iblocks)%Centroids(CenCount) = B1%Centroids(CenID) + BundlesMult(iblocks)%Centroids(CenCount)%CentID = [CBF_j, CBF_i] + BundlesMult(iblocks)%Centroids(CenCount)%TrajID = CenCount + write (fmiout, *) "We added Cen", CenID, "of B1 to block Bundle ", iblocks, "as Cen number", CenCount + + if (iblocks>1) then + !BlockCBFCount = 0 + !do jblocks = 1, iblocks-1 + ! BlockCBFCount = BlockCBFCount + NumCBFBlock(jblocks) + !end do + + ! still in the same block, so reuse BlockCBFCount from above + + BundlesMult(iblocks)%Centroids(CenCount)%CentID = [CBF_j - BlockCBFCount, & + CBF_i - BlockCBFCount] + write (fmiout, *) "Cen receives the number (must start at 1)" + end if + CenCount = CenCount + 1 + cycle + end if + + else + write (fmiout, *) "Index pair is not the one needed for the CentIDs in this block" + + end if + + end do + + write (fmiout, *) "___________________________" + + end do + + end do + + write (fmiout, *) "Do we have Centroids after copying them?", glzCentroids + do iblocks = 1, BundlesMultDim + if (BundlesMult(iblocks)%NCBFs > 1) then + do iTraj = 1, (((BundlesMult(iblocks)%NCBFs - 1) * BundlesMult(iblocks)%NCBFs) / 2) + write (fmiout, *) "Centroid number ", iTraj + write (fmiout, *) "is centroid to trajectories ", BundlesMult(iblocks)%Centroids(iTraj)%CentID + write (fmiout, *) "And has position: ", BundlesMult(iblocks)%Centroids(iTraj)%Particle(1)%get_pos() + end do + else + write (fmiout, *) "This block only has one CBF: ", iblocks, ". Therefore no centroids." + end if + end do + + write (fmiout, *) " ----------------------------------------------------" + + write (fmiout, *) "What we filled into the separated bundle: " + do iblocks = 1, BundlesMultDim + do iTraj = 1, BundlesMult(iblocks)%NumTraj + + write (fmiout, *) "Info for trajectory with ID ", BundlesMult(iblocks)%Trajectory(iTraj)%TrajID + write (fmiout, *) "Triplet?", BundlesMult(iblocks)%Trajectory(iTraj)%triplet + write (fmiout, *) "Ms", BundlesMult(iblocks)%Trajectory(iTraj)%Ms + write (fmiout, *) "CBF", BundlesMult(iblocks)%Trajectory(iTraj)%CBF + + end do + end do + write (fmiout, *) " ----------------------------------------------------" + + ! TODO: check that all sep bundles indeed have the right number of Centroids + + end subroutine getMultBundles + +! Copy over trajectories after the selection from BundlesMult +! to the original bundle. Note that we must copy over +! also the trajectories marked for death. + subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode) + integer(kind=DefInt), intent(in) :: BundlesMultDim + integer(kind=DefInt), dimension(BundlesMultDim), intent(in) :: StoSelMode + type(T_TrajectoryBundle), intent(inout) :: B1 + type(T_TrajectoryBundle), intent(in) :: BundlesMult(BundlesMultDim) + !type(T_TrajectoryBundle), pointer :: BSS_i + !type(T_Trajectory), pointer :: T_i + integer(kind=DefInt) :: i, j, iTraj + integer(kind=DefInt) :: TrajCount, CBFCount + !iSingTraj, iTripTraj, iCBF + !integer(kind=DefInt) :: TrajIDSing, CBFIDSing, TrajIDTrip, CBFIDTrip + + ! I dont actually think we need to copy back the temp separated bundles + ! All we need to do is to adjust gliForceKill to contain the TrajIDs for B1 + ! not those referring to the temp separated bundles + + do i = 1, BundlesMultDim + if (StoSelMode(i) == 1 .or. StoSelMode(i) == 3) then + write (fmiout, *) "Killing happened in block Bundle ", i + write (fmiout, *) "Trajs need to be copied back to B1 to keep track of who died" + + do iTraj = 1, BundlesMult(i)%NumTraj + write (fmiout, *) " ------------------------" + write (fmiout, *) "Copying back Traj ", iTraj + write (fmiout, *) "Ms", BundlesMult(i)%Trajectory(iTraj)%Ms + write (fmiout, *) "CBF", BundlesMult(i)%Trajectory(iTraj)%CBF + write (fmiout, *) "What is the amplitude?", BundlesMult(i)%Trajectory(iTraj)%Amplitude + write (fmiout, *) " ------------------------" + write (fmiout, *) "Is this traj still the same in B1?", B1%Trajectory(iTraj)%TrajID + write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms + write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF + write (fmiout, *) "What is the amplitude?", B1%Trajectory(iTraj)%Amplitude + write (fmiout, *) " ------------------------" + end do + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "Now, actually copy and check again: " + write (fmiout, *) " ----------------------------------------------------" + + ! Copy back trajectories, adjust temp bundle TrajID, CBF back to refer to B1 + write (fmiout, *) "Copying trajs for block bundle ", i + + TrajCount = 0 + CBFCount = 0 + if (i>1) then ! Problem here: we reduce BundlesMultDim + ! when merging Coupled, so BundlesMultDim + ! does not represent the total number of + ! blocks anymore + ! can we store the Coupled(B1) and use here? + do j = 1, i - 1 + TrajCount = TrajCount + BundlesMult(j)%NumTraj + CBFCount = CBFCount + BundlesMult(j)%NCBFs + end do + write (fmiout, *) "The", i-1, "previous bundles had: " + write (fmiout, *) TrajCount, "Trajs and", CBFcount, "CBFs" + + end if + + do iTraj = 1, BundlesMult(i)%NumTraj + if (i>1) then + B1%Trajectory(iTraj+TrajCount) = BundlesMult(i)%Trajectory(iTraj) !TODO: figure out what original TrajID, CBF ID was + B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF + BundlesMult(i-1)%NCBFs + end if + end do + + do iTraj = 1, BundlesMult(i)%NumTraj + write (fmiout, *) " ------------------------" + write (fmiout, *) "Copying back Traj ", iTraj + write (fmiout, *) "Ms", BundlesMult(i)%Trajectory(iTraj)%Ms + write (fmiout, *) "CBF", BundlesMult(i)%Trajectory(iTraj)%CBF + write (fmiout, *) "What is the amplitude?", BundlesMult(i)%Trajectory(iTraj)%Amplitude + write (fmiout, *) " ------------------------" + write (fmiout, *) "Is this traj still the same in B1?", B1%Trajectory(iTraj)%TrajID + write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms + write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF + write (fmiout, *) "What is the amplitude?", B1%Trajectory(iTraj)%Amplitude + write (fmiout, *) " ------------------------" + end do + end if + end do + + end subroutine copy_MultBundles_to_original_bundle + subroutine fill_state_bundles(B1, BundleSS) type(T_TrajectoryBundle), intent(in) :: B1 type(T_TrajectoryBundle), intent(inout) :: BundleSS(B1%NumStates) diff --git a/src/openfms.F90 b/src/openfms.F90 index 23e5271..e7bb90b 100644 --- a/src/openfms.F90 +++ b/src/openfms.F90 @@ -413,17 +413,22 @@ subroutine initialize_bundle(B1, NumTraj, NumParticles, NumStates, InitialState) DMTemp = 1.0 if (NumInitBasis > 1) then - T1(1)%Amplitude=DMTemp - !T1(1)%Amplitude=DMTemp / NumInitBasis ! for equally distributing initial amp + !T1(1)%Amplitude=DMTemp + T1(1)%Amplitude=DMTemp / NumInitBasis ! for equally distributing initial amp do ITraj = 2, NumTraj T1(ITraj) = T1(1) ! DH: This seems weird??? - T1(ITraj)%Amplitude=0.0 -! T1(ITraj)%Amplitude=DMTemp / NumInitBasis ! for equally distributing initial amp + !T1(ITraj)%Amplitude=0.0 + write (fmiout, *) "Theamplitude of traj ", ITraj, " will be set to", DMTemp / NumInitBasis + T1(ITraj)%Amplitude=DMTemp / NumInitBasis ! for equally distributing initial amp ! T1(ITraj)%Amplitude = DMTemp end do + do ITraj = 1, NumTraj + write (fmiout, *) "This is the amplitude of traj ", ITraj, " before adding to Bundle:", T1(ITraj)%Amplitude + end do + else T1(1)%Amplitude = DMTemp endif @@ -442,10 +447,14 @@ subroutine initialize_bundle(B1, NumTraj, NumParticles, NumStates, InitialState) call B1%add_traj(T1(iTraj)) end do - write (fmiout, *) "Trajectories in initially in Bundle: (after adding them)", B1%NumTraj - write (fmiout, *) "And their IDs: " - write (fmiout, *) "TrajID 1:", B1%Trajectory(1)%TrajID - write (fmiout, *) "TrajID 2:", B1%Trajectory(2)%TrajID + do ITraj = 1, NumTraj + write (fmiout, *) "This is the actual amplitude of traj ", ITraj, " :", B1%Trajectory(ITraj)%Amplitude + end do + + !write (fmiout, *) "Trajectories in initially in Bundle: (after adding them)", B1%NumTraj + !write (fmiout, *) "And their IDs: " + !write (fmiout, *) "TrajID 1:", B1%Trajectory(1)%TrajID + !write (fmiout, *) "TrajID 2:", B1%Trajectory(2)%TrajID do IParticle = 1, NumParticles call Particle(IParticle)%destroy() @@ -457,15 +466,15 @@ subroutine initialize_bundle(B1, NumTraj, NumParticles, NumStates, InitialState) end do deallocate (T1) - write (fmiout, *) "Trajectories in initially in Bundle: (after clean up step)", B1%NumTraj - write (fmiout, *) "And their IDs: " - write (fmiout, *) "TrajID 1:", B1%Trajectory(1)%TrajID - write (fmiout, *) "TrajID 2:", B1%Trajectory(2)%TrajID - write (fmiout, *) "And their Amplitudes: " - write (fmiout, *) "TrajID 1:", B1%Trajectory(1)%Amplitude - write (fmiout, *) "TrajID 2:", B1%Trajectory(2)%Amplitude - write (fmiout, *) "Norm of the Bundle: " - write (fmiout, *) FMS_Norm(B1) + !write (fmiout, *) "Trajectories in initially in Bundle: (after clean up step)", B1%NumTraj + !write (fmiout, *) "And their IDs: " + !write (fmiout, *) "TrajID 1:", B1%Trajectory(1)%TrajID + !write (fmiout, *) "TrajID 2:", B1%Trajectory(2)%TrajID + !write (fmiout, *) "And their Amplitudes: " + !write (fmiout, *) "TrajID 1:", B1%Trajectory(1)%Amplitude + !write (fmiout, *) "TrajID 2:", B1%Trajectory(2)%Amplitude + !write (fmiout, *) "Norm of the Bundle: " + !write (fmiout, *) FMS_Norm(B1) end subroutine initialize_bundle From 50d6fe72bd4ee4dc722a19b96d5f19d8210ec422 Mon Sep 17 00:00:00 2001 From: verab98 Date: Wed, 22 Apr 2026 13:04:09 +0100 Subject: [PATCH 06/10] Separated StoSel (still WiP) Death is now handled currectly (haha) TrajIDs of DeadTrajs are not skipped anymore when killing --- src/modules/BundleCalcsModule.f90 | 10 +- src/modules/OverlapModule.f90 | 6 +- src/modules/PropagationModule.f90 | 6 +- src/modules/SelectionModule.f90 | 166 ++++++++++++++++++++++++++++-- src/remove_dead.f90 | 1 + 5 files changed, 168 insertions(+), 21 deletions(-) diff --git a/src/modules/BundleCalcsModule.f90 b/src/modules/BundleCalcsModule.f90 index 367a9d3..10789b4 100644 --- a/src/modules/BundleCalcsModule.f90 +++ b/src/modules/BundleCalcsModule.f90 @@ -562,16 +562,16 @@ subroutine FMS_BuildHS(Bundle) ! TBF within the CBF. When the next Ms=2 TBF appears, we ! know that a new CBF is reached. if (glzCentroids .and. T_i%CBF /= T_j%CBF) then - write(fmiOut,*) "i, j: ", i,j - write(fmiOut,*) "IDi, IDj: ",T_i%StateID,T_j%StateID - write(fmiOut,*) "Msi, Msj: ", T_i%Ms,T_j%Ms - write(fmiOut,*) "CBFi, CBFj: ", T_i%CBF,T_j%CBF + !write(fmiOut,*) "i, j: ", i,j + !write(fmiOut,*) "IDi, IDj: ",T_i%StateID,T_j%StateID + !write(fmiOut,*) "Msi, Msj: ", T_i%Ms,T_j%Ms + !write(fmiOut,*) "CBFi, CBFj: ", T_i%CBF,T_j%CBF !! if( T_i%Ms.eq.2 .and. T_j%Ms.eq.2 ) then CBFi = T_i%CBF CBFj = T_j%CBF ICent = ((CBFi - 2) * (CBFi - 1)) / 2 + CBFj - write (fmiout, *) "This is what we are trying to use as CentID in BundleCalcs: ", ICent + !write (fmiout, *) "This is what we are trying to use as CentID in BundleCalcs: ", ICent T_ij => Bundle%Centroids(ICent) !write (fmiout, *) "Position", Bundle%Centroids(ICent)%Particle(1)%get_pos() !! endif diff --git a/src/modules/OverlapModule.f90 b/src/modules/OverlapModule.f90 index 51f6b20..4c682bd 100644 --- a/src/modules/OverlapModule.f90 +++ b/src/modules/OverlapModule.f90 @@ -626,9 +626,9 @@ function overlap_V_trajectory(T_i, T_j, T_c, S_ij_precalc) result(V) !bfec if (glzCentroids .and. (CBFi /= CBFj)) then ic = T_c%CentID(1); jc = T_c%CentID(2) - write (fmiout, *) " ------------------------------------------------------------ " - write (fmiout, *) "In the OverlapModule: getting ic and jc", ic, jc - write (fmiout, *) " ------------------------------------------------------------ " + !write (fmiout, *) " ------------------------------------------------------------ " + !write (fmiout, *) "In the OverlapModule: getting ic and jc", ic, jc + !write (fmiout, *) " ------------------------------------------------------------ " end if is = T_i%StateID; js = T_j%StateID diff --git a/src/modules/PropagationModule.f90 b/src/modules/PropagationModule.f90 index 289b551..9d19d40 100644 --- a/src/modules/PropagationModule.f90 +++ b/src/modules/PropagationModule.f90 @@ -3,7 +3,8 @@ module PropagationModule use FMSModule use GlobalModule, only: DefInt, DefReal, DefComp, DP5, FPZero, & glzConstrain, glzCentroids, glcIntegType, glzFullyCoupled, & - FMS_DieError + FMS_DieError, & + fmiOut use TrajectoryModule use TrajectoryCalcsModule, only: FMS_GetForce, FMS_PhaseDot, FMS_PosDot, FMS_MomDot use BundleModule @@ -69,7 +70,6 @@ subroutine FMS_PropSplitOperator(Bundle, TimeStep) end do ! Renormalize if multi-state problem - if (Bundle%Trajectory(1)%NumStates /= 1) then call FMS_Renormalize(Bundle, OldNorms, DNorm0) end if @@ -479,6 +479,8 @@ subroutine FMS_Renormalize(B1, OldNorms, DNorm0) real(kind=DefReal) :: DNorm integer(kind=DefInt) :: IState, ITraj + write (fmiout, *) "This is OldNorms: ", OldNorms + if (B1%NumStates /= LastDimension) then if (LastDimension /= 0) deallocate (NewNorms) allocate (NewNorms(B1%NumStates)) diff --git a/src/modules/SelectionModule.f90 b/src/modules/SelectionModule.f90 index 031d370..2e96def 100644 --- a/src/modules/SelectionModule.f90 +++ b/src/modules/SelectionModule.f90 @@ -6,6 +6,7 @@ module SelectionModule use BundleCalcsModule, only: FMS_bH, FMS_Norm use OverlapModule, only: overlap use RandomModule, only: fms_ranb + use PropagationModule, only: FMS_Renormalize implicit none private @@ -71,13 +72,13 @@ subroutine FMS_StochasticCollapse(B1) integer(kind=DefInt), dimension(B1%NumTraj + 1, B1%NumTraj + 1) :: blocktrajid integer(kind=DefInt), dimension(B1%NumTraj + 1, B1%NumTraj + 1) :: MergedBlockTrajID, MergedSepBlockTrajID integer(kind=DefInt), dimension(B1%NumTraj + 1) :: ntrajblock - integer(kind=DefInt) :: nblockMerged, DeadID + integer(kind=DefInt) :: nblockMerged, DeadID, NumDeadTrajCurr, NumDeadTrajSave integer(kind=DefInt), dimension(50) :: SaveForceKill ! sepSS: Shall we perform multiplet separated stochastic selection integer(kind=DefInt), allocatable :: StoSelMode(:), NumTrajBlock(:) logical :: isSing, isTrip, isSingTrip ! --- Added 4 sep SS end --- - integer(kind=DefInt) :: iTraj, i, j, nblock + integer(kind=DefInt) :: iTraj, i, j, k, nblock, totDeadCount ! AIMSWISS: Shall we perform stochastic selection logical :: performSelection ! AIMSWISS: Current selection time @@ -133,13 +134,14 @@ subroutine FMS_StochasticCollapse(B1) else if (NTrip /= 0) then + NumDeadTrajSave = B1%NumDeadTraj + ! (1) getting Coupled matrix to identify separated S only, or T only blocks ! - get the total converged coupled matrix, and the number of blocks write (fmiout, *) "-----------------------------------------------------------" write (fmiout, *) "Building coupled matrix for total Bundle (B1, unseparated)" Coupled = 0 - !write (fmiout, *) "This is the total, converged coupled matrix of size ", size(Coupled) call FMS_BuildCoupled(B1, Coupled, selectionTime) write (fmiout, *) "This is the total, converged coupled matrix of size ", size(Coupled) do iTraj = 1, B1%NumTraj @@ -223,11 +225,20 @@ subroutine FMS_StochasticCollapse(B1) end do end do + MergedBlockTrajID = blocktrajid nblockMerged = nblock + ! Generate 2nd version of MergedBlockTrajID w/ indices ref to what will become the sep bundles ! This is the MergedSepBlockTrajID matrx (sorry complicated.. :( ) + ! In MergedSepBlockTrajID every column restarts from 1, + ! while in MergedBlockTrajID indices continue + ! ------------------------------------------------------------------------------------------------ ! TODO: this doesn't seem robust at all, go back later and improve + ! edit: seems a little more robust now, at least (still problematic: + ! - is it fine to only adjust indices after the + ! first column = from the second block onwards??) + MergedSepBlockTrajID = MergedBlockTrajID do i = 1, B1%NumTraj + 1 do j = 2, B1%NumTraj + 1 @@ -236,18 +247,53 @@ subroutine FMS_StochasticCollapse(B1) end do end do + ! Make sure that indices in MergedBlockTrajID leave out indices of trajs that died before + totDeadCount = 0 + do i = 1, B1%NumTraj + 1 + do j = 2, B1%NumTraj + 1 + if (MergedBlockTrajID(i,j) == 0) cycle + do k = 1, B1%NumDeadTraj + if (MergedBlockTrajID(i,j) == B1%DeadTraj(k)%TrajID) then + write (fmiout, *) "This index belongs to a dead traj: ", MergedBlockTrajID(i,j), B1%DeadTraj(k)%TrajID + write (fmiout, *) "totDeadCount increased from ", totDeadCount + totDeadCount = totDeadCount + 1 + write (fmiout, *) "to... ", totDeadCount + end if + end do + MergedBlockTrajID(i,j) = MergedBlockTrajID(i,j) + totDeadCount + end do + end do + + + ! ----------------------------- !!! -------------------------------------------------- + ! 21/04/26: Dead Traj IDs need to be incorporated to MergedBlockTrajID + ! + ! 22/04/26: Check with unsep SS: does BlockTrajID account for Dead Trajs? + ! If so, why does it not in my case? + ! ----------------------------- !!! -------------------------------------------------- + write (fmiout, *) "This is the S, T, S/T StoSel mode merged blocktrajid: " do j = 1, B1%NumTraj + 1 write (fmiout, *) blocktrajid(j,:) end do write (fmiout, *) "And these are our StoSel (updated) modes: ", StoSelMode - write (fmiout, *) "This is the merged sepblocktrajid with what will be sep bundle IDs: " + write (fmiout, *) "This is the MergedBlockTrajID: " + do j = 1, B1%NumTraj + 1 + write (fmiout, *) MergedBlockTrajID(j,:) + end do + + write (fmiout, *) "This is the MergedSepBlockTrajID with what will be sep bundle IDs: " do j = 1, B1%NumTraj + 1 write (fmiout, *) MergedSepBlockTrajID(j,:) end do write (fmiout, *) "And these are our StoSel (updated) modes: ", StoSelMode + ! ----------------------------- !!! -------------------------------------------------- + + ! Now we have everything for figuring out in which of the separated bundles selection + ! needs to happen, and translating for back from the separated bundles to the total B1 + ! (2) Create and fill temp separated bundles ! Determine how many block bundles we need to create, based on merged blocktrajid: @@ -271,6 +317,7 @@ subroutine FMS_StochasticCollapse(B1) call getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) ! (3) Perform SS on the separated bundles + NumDeadTrajCurr = 0 do i = 1, nblock if (StoSelMode(i) == 1 .or. StoSelMode(i) == 3) then @@ -323,6 +370,18 @@ subroutine FMS_StochasticCollapse(B1) write (fmiout, *) gliForceKill end if + if (B1%NumDeadTraj > 0) then + + do iTraj = 1, B1%NumDeadTraj + write (fmiout, *) "IDs of dead traj ", iTraj, ": ", B1%DeadTraj(iTraj)%TrajID + end do + + !if (B1%NumDeadTraj >= 2) then + ! call FMS_DieError("That's enough for now...") + !end if + + end if + else write (fmiout, *) "No stochastic selection because block Bundle only has one CBF" end if @@ -337,6 +396,7 @@ subroutine FMS_StochasticCollapse(B1) end do + !write (fmiout, *) "Check that Centroids of B1 remain unaffected" !write (fmiout, *) "Part 1: B1 Centroids before" !do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) @@ -347,10 +407,20 @@ subroutine FMS_StochasticCollapse(B1) ! (4) Set appropriate amplitudes to zero in the original bundle ! If selection happened, copy information which trajs died - if (any(gliForceKill(:) /= 0)) then - call copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode) - write (fmiout, *) "States were copied back without error" - end if + + !if (any(gliForceKill(:) /= 0)) then ! TODO: use something else here + ! otherwise copy_MultBundles_to_original_bundle always gets called as soon as + ! we have a number in gliForceKill + ! How to determine if kill happened? + ! --- + ! For now, go with setting gliForceKill back to zero at the start of each + ! StoSel (but not sure if that breaks things somewhere...?) + + + call copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode) + write (fmiout, *) "States were copied back without error" + + !end if !write (fmiout, *) "Check that Centroids of B1 remain unaffected" !write (fmiout, *) "Part 2: B1 Centroids after" @@ -1402,7 +1472,11 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) end do write (fmiout, *) " ----------------------------------------------------" - ! TODO: check that all sep bundles indeed have the right number of Centroids + ! Normalise the separate bundles before SS (and keep the og norm, amplitudes?) + + !do iblocks = 1, BundlesMultDim + ! call FMS_Renormalize(BundlesMult(i), + !end do end subroutine getMultBundles @@ -1417,7 +1491,9 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, !type(T_TrajectoryBundle), pointer :: BSS_i !type(T_Trajectory), pointer :: T_i integer(kind=DefInt) :: i, j, iTraj - integer(kind=DefInt) :: TrajCount, CBFCount + integer(kind=DefInt) :: TrajCount, CBFCount, DeadCount + integer(kind=DefInt), dimension(50) :: DeadTrajIDs + logical :: hitDead !iSingTraj, iTripTraj, iCBF !integer(kind=DefInt) :: TrajIDSing, CBFIDSing, TrajIDTrip, CBFIDTrip @@ -1425,6 +1501,50 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, ! All we need to do is to adjust gliForceKill to contain the TrajIDs for B1 ! not those referring to the temp separated bundles + do i = 1, B1%NumTraj + write (fmiout, *) "These trajs are in the original B1 (before copying back)" + write (fmiout, *) "TrajID: ", B1%Trajectory(i)%TrajID, "DeadTime: ", B1%Trajectory(i)%DeadTime + write (fmiout, *) "IsDead? (checking if DeadTime < CurrentTime)", B1%Trajectory(i)%DeadTime < B1%CurrentTime + end do + + ! How to get the information on which Trajs are dead? + !write (fmiout, *) "These trajs are dead in the original B1 (before copying back)" + !do i = 1, B1%NumDeadTraj + ! write (fmiout, *) "TrajID: ", B1%DeadTraj(i)%TrajID, "DeadTime: ", B1%DeadTraj(i)%DeadTime + !end do + + ! Try a different approach to what is (or will be) commented out further down + + DeadCount = 0 + DeadTrajIDs = 0 + do iTraj = 1, B1%NumTraj ! Does this loop always account for all dead and alive trajs??? + ! Even after several kills in the Bundle? + ! Do we need to add B1%NumDeadTraj somehow asp? + write (fmiout, *) " ------------------------" + write (fmiout, *) "Current iTraj: ", iTraj + + if (any(B1%Trajectory(:)%TrajID == iTraj)) then + write(fmiout, *) "This traj is alive" + end if + + if (any(B1%DeadTraj(:)%TrajID == iTraj)) then + write(fmiout, *) "This traj is dead" + ! This TrajID needs to be remembered and 'left blank' when temp separated Bundles are copied back + DeadCount = DeadCount + 1 + DeadTrajIDs(DeadCount) = iTraj + end if + + if (any(B1%DeadTraj(:)%TrajID == iTraj) .and. any(B1%Trajectory(:)%TrajID == iTraj)) then + write(fmiout, *) "This traj is dead and alive" + call FMS_DieError("This shouldn't happen") + end if + + write (fmiout, *) " ------------------------" + end do + write (fmiout, *) "These are our dead trajectories: ", DeadTrajIDs + + hitDead = .false. + do i = 1, BundlesMultDim if (StoSelMode(i) == 1 .or. StoSelMode(i) == 3) then write (fmiout, *) "Killing happened in block Bundle ", i @@ -1469,8 +1589,23 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, do iTraj = 1, BundlesMult(i)%NumTraj if (i>1) then B1%Trajectory(iTraj+TrajCount) = BundlesMult(i)%Trajectory(iTraj) !TODO: figure out what original TrajID, CBF ID was - B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + + if (any(BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj == DeadTrajIDs(:))) then + write (fmiout, *) "We hit the dead traj (iTraj)", & + BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + write (fmiout, *) "Adding the DeadCount ", DeadCount, "from now on" + hitDead = .true. + end if + + if (hitDead) then + B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + & + DeadCount + else + B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + end if B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF + BundlesMult(i-1)%NCBFs + write (fmiout, *) "Traj of index ", iTraj, "was copied over and received ID, ", & + B1%Trajectory(iTraj+TrajCount)%TrajID end if end do @@ -1490,6 +1625,15 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, end if end do + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "These trajs are in B1 (after copying back)" + do i = 1, B1%NumTraj + write (fmiout, *) "TrajID: ", B1%Trajectory(i)%TrajID, "CBF: ", B1%Trajectory(i)%CBF + write (fmiout, *) "Ms: ", B1%Trajectory(i)%Ms + write (fmiout, *) "Amp: ", B1%Trajectory(i)%Amplitude + end do + write (fmiout, *) " ----------------------------------------------------" + end subroutine copy_MultBundles_to_original_bundle subroutine fill_state_bundles(B1, BundleSS) diff --git a/src/remove_dead.f90 b/src/remove_dead.f90 index f7144ba..70574ec 100644 --- a/src/remove_dead.f90 +++ b/src/remove_dead.f90 @@ -125,6 +125,7 @@ subroutine FMS_RemoveDead(B1) ndead = ndead + 1 if (ForceDead) then write (fmiOut, '(a,i0,a,i0)') "** Force Killing trajectory ", TrajID, " on state ", StateID + write (fmiout, *) "Based on gliForceKill: ", gliForceKill elseif (OnIgnoreState) then write (fmiOut, '(a,i0,a,i0)') "** Killing trajectory ", TrajID, " on state ", StateID else From caecbe7b009da0ed3d30bd8b3133fc4b6c5ca71d Mon Sep 17 00:00:00 2001 From: verab98 Date: Wed, 29 Apr 2026 16:38:30 +0100 Subject: [PATCH 07/10] What looks like a version that is working - Amplitudes are normalised to separated block and renormalised to total bundle after selection - Issue of no normalisation if selection happens in first separated block fixed - Index of trajectories to kill only changed if selection happened (added IsSelection to perform_stochastic_selection for this) - Commented out a bunch of print statements --- src/dynamics.f90 | 46 +- src/modules/PropagationModule.f90 | 2 +- src/modules/SelectionModule.f90 | 749 +++++++++++++++++++----------- src/remove_dead.f90 | 24 +- 4 files changed, 513 insertions(+), 308 deletions(-) diff --git a/src/dynamics.f90 b/src/dynamics.f90 index b6a038d..b99272f 100644 --- a/src/dynamics.f90 +++ b/src/dynamics.f90 @@ -80,32 +80,32 @@ recursive subroutine FMS_Dynamics(Bundle, Timestep, GoalTime) !Stochastic collapse nuclear wavefunction !(RemoveDead must be called immediately after) call FMS_StochasticCollapse(Bundle) - write (fmiout, *) "-----------------------------------------------------------" - write (fmiout, *) "FMS_StochasticCollapse was successfully called, leaving now" - write (fmiout, *) "Now the time is: ", Bundle%CurrentTime - write (fmiout, *) "And these are our trajectories: (NumTraj, NCBFs)", Bundle%NumTraj, Bundle%NCBFs - write (fmiout, *) "And these are their centroids: " - do iTraj = 1, (((Bundle%NCBFs - 1) * Bundle%NCBFs) / 2) - write (fmiout, *) iTraj - write (fmiout, *) "is Cen to trajectories: ", Bundle%Centroids(iTraj)%CentID - write (fmiout, *) "Position", Bundle%Centroids(iTraj)%Particle(1)%get_pos() - end do - write (fmiout, *) "-----------------------------------------------------------" + !write (fmiout, *) "-----------------------------------------------------------" + !write (fmiout, *) "FMS_StochasticCollapse was successfully called, leaving now" + !write (fmiout, *) "Now the time is: ", Bundle%CurrentTime + !write (fmiout, *) "And these are our trajectories: (NumTraj, NCBFs)", Bundle%NumTraj, Bundle%NCBFs + !write (fmiout, *) "And these are their centroids: " + !do iTraj = 1, (((Bundle%NCBFs - 1) * Bundle%NCBFs) / 2) + ! write (fmiout, *) iTraj + ! write (fmiout, *) "is Cen to trajectories: ", Bundle%Centroids(iTraj)%CentID + ! write (fmiout, *) "Position", Bundle%Centroids(iTraj)%Particle(1)%get_pos() + !end do + !write (fmiout, *) "-----------------------------------------------------------" end if - write (fmiout, *) "Now calling RemoveDead... " + !write (fmiout, *) "Now calling RemoveDead... " call FMS_RemoveDead(Bundle) if (Bundle%NumTraj == 0) return - write (fmiout, *) "-----------------------------------------------------------" - write (fmiout, *) "Did deleting dead parts of Bundle mess up anything?" - write (fmiout, *) "Now the time is: ", Bundle%CurrentTime - write (fmiout, *) "And these are our trajectories: (NumTraj, NCBFs)", Bundle%NumTraj, Bundle%NCBFs - write (fmiout, *) "And these are their centroids: " - do iTraj = 1, (((Bundle%NCBFs - 1) * Bundle%NCBFs) / 2) - write (fmiout, *) iTraj - write (fmiout, *) "is Cen to trajectories: ", Bundle%Centroids(iTraj)%CentID - write (fmiout, *) "Position", Bundle%Centroids(iTraj)%Particle(1)%get_pos() - end do - write (fmiout, *) "-----------------------------------------------------------" + !write (fmiout, *) "-----------------------------------------------------------" + !write (fmiout, *) "Did deleting dead parts of Bundle mess up anything?" + !write (fmiout, *) "Now the time is: ", Bundle%CurrentTime + !write (fmiout, *) "And these are our trajectories: (NumTraj, NCBFs)", Bundle%NumTraj, Bundle%NCBFs + !write (fmiout, *) "And these are their centroids: " + !do iTraj = 1, (((Bundle%NCBFs - 1) * Bundle%NCBFs) / 2) + ! write (fmiout, *) iTraj + ! write (fmiout, *) "is Cen to trajectories: ", Bundle%Centroids(iTraj)%CentID + ! write (fmiout, *) "Position", Bundle%Centroids(iTraj)%Particle(1)%get_pos() + !end do + !write (fmiout, *) "-----------------------------------------------------------" end if ! Write output here only if user requested output at every step diff --git a/src/modules/PropagationModule.f90 b/src/modules/PropagationModule.f90 index 9d19d40..0dbf3d5 100644 --- a/src/modules/PropagationModule.f90 +++ b/src/modules/PropagationModule.f90 @@ -138,7 +138,7 @@ subroutine FMS_PropQCVV(Bundle, TimeStep) id = Bundle%Trajectory(j)%TrajID amp = Bundle%Trajectory(j)%Amplitude ms = Bundle%Trajectory(j)%Ms - ! som=Bundle%Trajectory(j)%ElecStruc%SOMat + ! som=Bundle%Trajectory(j)%ElecStruc%SOMat ! Why is SOMat not copied??? cohist = Bundle%Trajectory(j)%CoupHist call Bundle%Trajectory(j)%copy_from(Bundle%Trajectory(i)) diff --git a/src/modules/SelectionModule.f90 b/src/modules/SelectionModule.f90 index 2e96def..b94c01d 100644 --- a/src/modules/SelectionModule.f90 +++ b/src/modules/SelectionModule.f90 @@ -3,10 +3,9 @@ module SelectionModule use GlobalModule use BundleModule use TrajectoryModule - use BundleCalcsModule, only: FMS_bH, FMS_Norm + use BundleCalcsModule, only: FMS_bH, FMS_Norm, FMS_Branching use OverlapModule, only: overlap use RandomModule, only: fms_ranb - use PropagationModule, only: FMS_Renormalize implicit none private @@ -74,11 +73,14 @@ subroutine FMS_StochasticCollapse(B1) integer(kind=DefInt), dimension(B1%NumTraj + 1) :: ntrajblock integer(kind=DefInt) :: nblockMerged, DeadID, NumDeadTrajCurr, NumDeadTrajSave integer(kind=DefInt), dimension(50) :: SaveForceKill + real(kind=DefReal), allocatable :: SaveNorms(:) + real(kind=DefReal), allocatable :: SaveTrajAmps(:,:) ! sepSS: Shall we perform multiplet separated stochastic selection integer(kind=DefInt), allocatable :: StoSelMode(:), NumTrajBlock(:) logical :: isSing, isTrip, isSingTrip + logical :: IsSelection ! --- Added 4 sep SS end --- - integer(kind=DefInt) :: iTraj, i, j, k, nblock, totDeadCount + integer(kind=DefInt) :: iTraj, i, j, k, nblock_outer, totDeadCount ! AIMSWISS: Shall we perform stochastic selection logical :: performSelection ! AIMSWISS: Current selection time @@ -136,25 +138,28 @@ subroutine FMS_StochasticCollapse(B1) NumDeadTrajSave = B1%NumDeadTraj + !write (fmiout, *) "Asking this question once and for all: do we call FMS_StochasticCollapse at each time step??" + ! yes we do... + ! (1) getting Coupled matrix to identify separated S only, or T only blocks ! - get the total converged coupled matrix, and the number of blocks - write (fmiout, *) "-----------------------------------------------------------" - write (fmiout, *) "Building coupled matrix for total Bundle (B1, unseparated)" + !write (fmiout, *) "-----------------------------------------------------------" + !write (fmiout, *) "Building coupled matrix for total Bundle (B1, unseparated)" Coupled = 0 call FMS_BuildCoupled(B1, Coupled, selectionTime) - write (fmiout, *) "This is the total, converged coupled matrix of size ", size(Coupled) - do iTraj = 1, B1%NumTraj - write (fmiout, *) Coupled(iTraj,:) - end do + !write (fmiout, *) "This is the total, converged coupled matrix of size ", size(Coupled) + !do iTraj = 1, B1%NumTraj + ! write (fmiout, *) Coupled(iTraj,:) + !end do - write (fmiout, *) "Get number of blocks for total Bundle (B1, unseparated)" + !write (fmiout, *) "Get number of blocks for total Bundle (B1, unseparated)" blocktrajid = 0 ntrajblock = 0 - nblock = 0 - call FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, ntrajblock, nblock) - write (fmiout, *) "Number of Blocks: ", nblock - write (fmiout, *) "-----------------------------------------------------------" + nblock_outer = 0 + call FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, ntrajblock, nblock_outer) + !write (fmiout, *) "Number of Blocks: ", nblock_outer + !write (fmiout, *) "-----------------------------------------------------------" ! - merge the total coupled matrix according to ! - Singlet only --> perform StoSel @@ -164,15 +169,15 @@ subroutine FMS_StochasticCollapse(B1) isTrip = .false. isSingTrip = .false. - write (fmiout, *) "-----------------------------------------------------------" - write (fmiout, *) "Deciding StoSel mode based on total Bundle" + !write (fmiout, *) "-----------------------------------------------------------" + !write (fmiout, *) "Deciding StoSel mode based on total Bundle" - allocate (StoSelMode(nblock)) - allocate (NumTrajBlock(nblock)) + allocate (StoSelMode(nblock_outer)) + allocate (NumTrajBlock(nblock_outer)) ! Get isSIng, isTrip, isSingTrip for each block in the original, total Coupled matrix NumTrajBlock = 0 - do i = 1, nblock + do i = 1, nblock_outer isSing = .false. isTrip = .false. isSingTrip = .false. @@ -198,19 +203,19 @@ subroutine FMS_StochasticCollapse(B1) call FMS_DieError("No valid Selection mode could be determined") end if - write (fmiout, *) "For block ", i, "StoSel mode is: " - write (fmiout, *) "isSing is ", isSing, ", isTrip is ", isTrip, ", and... isSingTrip is", isSingTrip - write (fmiout, *) "This is for the trajs ", blocktrajid(:,i) - write (fmiout, *) "This is the current StoSel mode: ", StoSelMode(i) + !write (fmiout, *) "For block ", i, "StoSel mode is: " + !write (fmiout, *) "isSing is ", isSing, ", isTrip is ", isTrip, ", and... isSingTrip is", isSingTrip + !write (fmiout, *) "This is for the trajs ", blocktrajid(:,i) + !write (fmiout, *) "This is the current StoSel mode: ", StoSelMode(i) end do ! Merge blocks depending on isSIng, isTrip ! (No action for isSingTrip because not to be touched in StoSel - do i = 1, nblock + do i = 1, nblock_outer if (StoSelMode(i) == 0) cycle - do j = i+1, nblock + do j = i+1, nblock_outer if (StoSelMode(i) == StoSelMode(j)) then ! Copy over block j TrajIDs to i column of blocktrajid @@ -227,7 +232,10 @@ subroutine FMS_StochasticCollapse(B1) end do MergedBlockTrajID = blocktrajid - nblockMerged = nblock + nblockMerged = nblock_outer ! TODO: For now, nblockMerged is not used, but do we need to account for change + ! in number of blocks after merging (so, nblockMerged < nblock_outer depending + ! on how many blocks we merged??? + ! This is actually done two steps later I think ! Generate 2nd version of MergedBlockTrajID w/ indices ref to what will become the sep bundles ! This is the MergedSepBlockTrajID matrx (sorry complicated.. :( ) @@ -254,10 +262,10 @@ subroutine FMS_StochasticCollapse(B1) if (MergedBlockTrajID(i,j) == 0) cycle do k = 1, B1%NumDeadTraj if (MergedBlockTrajID(i,j) == B1%DeadTraj(k)%TrajID) then - write (fmiout, *) "This index belongs to a dead traj: ", MergedBlockTrajID(i,j), B1%DeadTraj(k)%TrajID - write (fmiout, *) "totDeadCount increased from ", totDeadCount + !write (fmiout, *) "This index belongs to a dead traj: ", MergedBlockTrajID(i,j), B1%DeadTraj(k)%TrajID + !write (fmiout, *) "totDeadCount increased from ", totDeadCount totDeadCount = totDeadCount + 1 - write (fmiout, *) "to... ", totDeadCount + !write (fmiout, *) "to... ", totDeadCount end if end do MergedBlockTrajID(i,j) = MergedBlockTrajID(i,j) + totDeadCount @@ -272,22 +280,22 @@ subroutine FMS_StochasticCollapse(B1) ! If so, why does it not in my case? ! ----------------------------- !!! -------------------------------------------------- - write (fmiout, *) "This is the S, T, S/T StoSel mode merged blocktrajid: " - do j = 1, B1%NumTraj + 1 - write (fmiout, *) blocktrajid(j,:) - end do - write (fmiout, *) "And these are our StoSel (updated) modes: ", StoSelMode + !write (fmiout, *) "This is the S, T, S/T StoSel mode merged blocktrajid: " + !do j = 1, B1%NumTraj + 1 + ! write (fmiout, *) blocktrajid(j,:) + !end do + !write (fmiout, *) "And these are our StoSel (updated) modes: ", StoSelMode - write (fmiout, *) "This is the MergedBlockTrajID: " - do j = 1, B1%NumTraj + 1 - write (fmiout, *) MergedBlockTrajID(j,:) - end do + !write (fmiout, *) "This is the MergedBlockTrajID: " + !do j = 1, B1%NumTraj + 1 + ! write (fmiout, *) MergedBlockTrajID(j,:) + !end do - write (fmiout, *) "This is the MergedSepBlockTrajID with what will be sep bundle IDs: " - do j = 1, B1%NumTraj + 1 - write (fmiout, *) MergedSepBlockTrajID(j,:) - end do - write (fmiout, *) "And these are our StoSel (updated) modes: ", StoSelMode + !write (fmiout, *) "This is the MergedSepBlockTrajID with what will be sep bundle IDs: " + !do j = 1, B1%NumTraj + 1 + ! write (fmiout, *) MergedSepBlockTrajID(j,:) + !end do + !write (fmiout, *) "And these are our StoSel (updated) modes: ", StoSelMode ! ----------------------------- !!! -------------------------------------------------- @@ -299,16 +307,16 @@ subroutine FMS_StochasticCollapse(B1) ! Determine how many block bundles we need to create, based on merged blocktrajid: ! and based on whether we need to carry out SS in the merged block or not BundlesMultDim = 0 - do i = 1, nblock + do i = 1, nblock_outer if (StoSelMode(i) /= 0) then BundlesMultDim = BundlesMultDim + 1 end if end do - nblock = BundlesMultDim - write (fmiout, *) "BundlesMultDim based on merged blocktrajid: ", nblock - write (fmiout, *) "nblock based on merged blocktrajid: ", nblock + nblockMerged = BundlesMultDim + !write (fmiout, *) "BundlesMultDim based on merged blocktrajid: ", BundlesMultDim + !write (fmiout, *) "nblock based on merged blocktrajid: ", nblockMerged - write (fmiout, *) "-----------------------------------------------------------" + !write (fmiout, *) "-----------------------------------------------------------" allocate (BundlesMult(BundlesMultDim)) @@ -316,9 +324,50 @@ subroutine FMS_StochasticCollapse(B1) ! with corresponding trajs from B1 call getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) + ! Set the norm of each temporary bundle to one (for SS) + ! ---> by just normalizing + ! and remember actual norm within total bundle (for continuing propagation after SS) + ! ---> SaveNorms + + write (fmiout, *) "What's the norm of the newly created bundles?" + allocate (SaveNorms(BundlesMultDim)) + allocate (SaveTrajAmps(BundlesMultDim, B1%NumStates)) + do i = 1, BundlesMultDim + write (fmiout, *) "Block ", i, "Norm: ", FMS_Norm(BundlesMult(i)) + SaveNorms(i) = FMS_Norm(BundlesMult(i)) + call FMS_Branching(BundlesMult(i), SaveTrajAmps(i,:)) + end do + write (fmiout, *) "Did we copy over to SaveNorms correctly?", SaveNorms + write (fmiout, *) "Did we copy over to SaveTrajAmps correctly?" + !do i = 1, BundlesMultDim + ! write(fmiout, *) SaveTrajAmps(i, :) + !end do + + write (fmiout, *) "-----------------------------------------------------------" + write (fmiout, *) "The B1 norm is: " + write (fmiout, *) "Norm: ", FMS_Norm(B1) + write (fmiout, *) "and the corresponding amplitudes are" + do iTraj = 1, B1%NumTraj + write (fmiout, *) "Trajectory ", iTraj, "AMplitude: ", B1%Trajectory(iTraj)%Amplitude + end do + + write (fmiout, *) "----------------------------" + write (fmiout, *) "Normalizing..." + do i = 1, BundlesMultDim + do iTraj = 1, BundlesMult(i)%NumTraj + BundlesMult(i)%Trajectory(iTraj)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude / sqrt(SaveNorms(i)) + end do + end do + write (fmiout, *) "After Normalizing prior to SS, the sep bundle norms are: " + do i = 1, BundlesMultDim + write (fmiout, *) "Block ", i, "Norm: ", FMS_Norm(BundlesMult(i)) + end do + write (fmiout, *) "----------------------------" + + ! (3) Perform SS on the separated bundles NumDeadTrajCurr = 0 - do i = 1, nblock + do i = 1, nblockMerged if (StoSelMode(i) == 1 .or. StoSelMode(i) == 3) then @@ -328,25 +377,24 @@ subroutine FMS_StochasticCollapse(B1) write (fmiout, *) "Performing stochastic selection for block Bundle", i, & "with", BundlesMult(i)%NCBFs, "CBFs" - call perform_stochastic_selection(BundlesMult(i), selectionTime) ! TODO: make sure the right Trajs - ! get killed (currently: gliForceKill - ! uses TrajID based on sep bundle 2 - ! kill in B1 + IsSelection = .false. + call perform_stochastic_selection(BundlesMult(i), selectionTime, IsSelection) - if (i>1) then - write (fmiout, *) "Adjusting gliForceKill using MergedSepBlockTrajID (this is the full matrix): " - do j = 1, B1%NumTraj + 1 - write (fmiout, *) MergedSepBlockTrajID(j,:) - end do + ! Adjust gliForceKill only if selection actually happened + if (i>1 .and. IsSelection) then + !write (fmiout, *) "Adjusting gliForceKill using MergedSepBlockTrajID (this is the full matrix): " + !do j = 1, B1%NumTraj + 1 + ! write (fmiout, *) MergedSepBlockTrajID(j,:) + !end do - write (fmiout, *) "and MergedBlockTrajID (this is the full matrix): " - do j = 1, B1%NumTraj + 1 - write (fmiout, *) MergedBlockTrajID(j,:) - end do + !write (fmiout, *) "and MergedBlockTrajID (this is the full matrix): " + !do j = 1, B1%NumTraj + 1 + ! write (fmiout, *) MergedBlockTrajID(j,:) + !end do write (fmiout, *) "gliForcKill for block Bundle", i write (fmiout, *) gliForceKill - write (fmiout, *) "adjusting to B1 indiced..." + !write (fmiout, *) "adjusting to B1 indiced..." SaveForceKill = gliForceKill @@ -354,12 +402,12 @@ subroutine FMS_StochasticCollapse(B1) do j = 1, B1%NumTraj + 1 if (MergedSepBlockTrajID(iTraj,j) == 0) cycle do DeadID = 1, size(SaveForceKill) - write (fmiout, *) "many numbers (?)", SaveForceKill(DeadID) + !write (fmiout, *) "many numbers (?)", SaveForceKill(DeadID) if (MergedSepBlockTrajID(iTraj,j) == SaveForceKill(DeadID)) then - write (fmiout, *) "Match for gliForceKill: ", SaveForceKill(DeadID), "(this should be 4)" - write (fmiout, *) "with MergedSepblock: ", MergedSepBlockTrajID(iTraj,j) - write (fmiout, *) "i and j are: ", iTraj, j, "(should be 2, 4)" - write (fmiout, *) "Corresponding MergedBlock: ", MergedBlockTrajID(iTraj,j), "(this should be 6)" + !write (fmiout, *) "Match for gliForceKill: ", SaveForceKill(DeadID), "(this should be 4)" + !write (fmiout, *) "with MergedSepblock: ", MergedSepBlockTrajID(iTraj,j) + !write (fmiout, *) "i and j are: ", iTraj, j, "(should be 2, 4)" + !write (fmiout, *) "Corresponding MergedBlock: ", MergedBlockTrajID(iTraj,j), "(this should be 6)" gliForceKill(DeadID) = MergedBlockTrajID(iTraj,j) end if end do @@ -370,17 +418,17 @@ subroutine FMS_StochasticCollapse(B1) write (fmiout, *) gliForceKill end if - if (B1%NumDeadTraj > 0) then + !if (B1%NumDeadTraj > 0) then - do iTraj = 1, B1%NumDeadTraj - write (fmiout, *) "IDs of dead traj ", iTraj, ": ", B1%DeadTraj(iTraj)%TrajID - end do + ! do iTraj = 1, B1%NumDeadTraj + ! write (fmiout, *) "IDs of dead traj ", iTraj, ": ", B1%DeadTraj(iTraj)%TrajID + ! end do - !if (B1%NumDeadTraj >= 2) then - ! call FMS_DieError("That's enough for now...") - !end if + ! !if (B1%NumDeadTraj >= 2) then + ! ! call FMS_DieError("That's enough for now...") + ! !end if - end if + !end if else write (fmiout, *) "No stochastic selection because block Bundle only has one CBF" @@ -397,13 +445,42 @@ subroutine FMS_StochasticCollapse(B1) end do - !write (fmiout, *) "Check that Centroids of B1 remain unaffected" - !write (fmiout, *) "Part 1: B1 Centroids before" - !do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) - ! write (fmiout, *) "Centroid number", iTraj - ! write (fmiout, *) "Is centroid to trajectories", B1%Centroids(iTraj)%CentID - ! write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() - !end do + write (fmiout, *) "-----------------------------------------------------------" + write (fmiout, *) "Check that Centroids of B1 remain unaffected" + write (fmiout, *) "Part 1: B1 Centroids before" + do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) + write (fmiout, *) "Centroid number", iTraj + write (fmiout, *) "Is centroid to trajectories", B1%Centroids(iTraj)%CentID + write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() + end do + write (fmiout, *) "-----------------------------------------------------------" + + ! Set the norm of each temporary bundle back to SaveNorms + ! and renormalize adjusting traj amplitudes + ! ---> call FMS_Renormalize + write (fmiout, *) "----------------------------" + write (fmiout, *) "After Selection, the sep bundle norms are: " + do i = 1, BundlesMultDim + write (fmiout, *) "Block ", i, "Norm: ", FMS_Norm(BundlesMult(i)) + end do + + write (fmiout, *) "Renormalizing..." + do i = 1, BundlesMultDim + call FMS_Renormalize(BundlesMult(i), SaveTrajAmps(i,:), SaveNorms(i)) + end do + write (fmiout, *) "----------------------------" + write (fmiout, *) "After Renormalizing, the sep bundle norms are: " + do i = 1, BundlesMultDim + write (fmiout, *) "Block ", i, "Norm: ", FMS_Norm(BundlesMult(i)) + end do + write (fmiout, *) "and the corresponding amplitudes are" + do i = 1, BundlesMultDim + write (fmiout, *) "Block: ", i + do iTraj = 1, BundlesMult(i)%NumTraj + write (fmiout, *) "Trajectory ", iTraj, "AMplitude: ", BundlesMult(i)%Trajectory(iTraj)%Amplitude + end do + end do + write (fmiout, *) "----------------------------" ! (4) Set appropriate amplitudes to zero in the original bundle ! If selection happened, copy information which trajs died @@ -418,10 +495,19 @@ subroutine FMS_StochasticCollapse(B1) call copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode) - write (fmiout, *) "States were copied back without error" + !write (fmiout, *) "States were copied back without error" - !end if + write (fmiout, *) "After copy_MultBundles_to_original_Bundle, the B1 norm is: " + write (fmiout, *) "Norm: ", FMS_Norm(B1) + write (fmiout, *) "and the corresponding amplitudes are" + do iTraj = 1, B1%NumTraj + write (fmiout, *) "Trajectory ", iTraj, "AMplitude: ", B1%Trajectory(iTraj)%Amplitude + end do + write (fmiout, *) "-----------------------------------------------------------" + !end if + + !write (fmiout, *) " ----------------------------------------------------" !write (fmiout, *) "Check that Centroids of B1 remain unaffected" !write (fmiout, *) "Part 2: B1 Centroids after" !do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) @@ -429,11 +515,13 @@ subroutine FMS_StochasticCollapse(B1) ! write (fmiout, *) "Is centroid to trajectories", B1%Centroids(iTraj)%CentID ! write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() !end do + !write (fmiout, *) " ----------------------------------------------------" ! Destroy the temporary block bundles - do i = 1, nblock + do i = 1, nblockMerged call BundlesMult(i)%destroy() end do + deallocate(SaveNorms, SaveTrajAmps, BundlesMult) !call FMS_DieError("Now, please go back and reconsider your life choices!") @@ -458,7 +546,7 @@ end subroutine FMS_StochasticCollapse !! This is done by renormalizing population in collapsed block to one !! and then marking all other trajectories to be killed. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine perform_stochastic_selection(B1, selectionTime) + subroutine perform_stochastic_selection(B1, selectionTime, IsSelection) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - type(T_TrajectoryBundle), intent(inout) :: B1 ! Coupling matrix for TBF basis @@ -477,14 +565,16 @@ subroutine perform_stochastic_selection(B1, selectionTime) integer(kind=DefInt) :: iblockslct ! AIMSWISS: Current selection time real(kind=DefReal), intent(in) :: selectionTime +! sep SS: Did selection happen? + logical, optional, intent(inout) :: IsSelection - write (fmiout, *) "Number of traj used for dimension of Coupled: ", B1%NumTraj - write (fmiout, *) "meaning Coupled dimension is: ", "( ", B1%NumTraj, ", ", B1%NumTraj, " )" + !write (fmiout, *) "Number of traj used for dimension of Coupled: ", B1%NumTraj + !write (fmiout, *) "meaning Coupled dimension is: ", "( ", B1%NumTraj, ", ", B1%NumTraj, " )" - write (fmiout, *) "CentIDs in perform_stochastic_selection (iTraj, iCBF:" - do iTraj = 1, B1%NumTraj - write (fmiout, *) B1%Trajectory(iTraj)%TrajID, B1%Trajectory(iTraj)%CBF - end do + !write (fmiout, *) "CentIDs in perform_stochastic_selection (iTraj, iCBF:" + !do iTraj = 1, B1%NumTraj + ! write (fmiout, *) B1%Trajectory(iTraj)%TrajID, B1%Trajectory(iTraj)%CBF + !end do ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ! First, decompose the hamiltonian into a block diagonal representation ! work out the coupling matrix @@ -505,9 +595,13 @@ subroutine perform_stochastic_selection(B1, selectionTime) call FMS_DieError("ERROR in FMS_StochasticCollapse") end if ! If just a single block then no more work to be done + if (nblock == 1) write (fmiout, *) "Only one block. No selection!" if (nblock == 1) return ! At this point, more than one uncoupled block of trajectories, ! so we want to stochastically collapse to one of them + if (present(IsSelection)) then + IsSelection = .true. + end if ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ! Compute populations of each block @@ -687,9 +781,9 @@ subroutine FMS_BuildCoupled(B1, Coupled, selectionTime) real(kind=DefReal), intent(in) :: selectionTime integer(kind=DefInt) :: ntraj, i, j - write (fmiout, *) "---------------------------------------------------------------------------" - write (fmiout, *) " Inside FMS_BuildCoupled " - write (fmiout, *) "The number of traj we are using is: ", B1%NumTraj + !write (fmiout, *) "---------------------------------------------------------------------------" + !write (fmiout, *) " Inside FMS_BuildCoupled " + !write (fmiout, *) "The number of traj we are using is: ", B1%NumTraj ntraj = B1%NumTraj @@ -700,8 +794,8 @@ subroutine FMS_BuildCoupled(B1, Coupled, selectionTime) if (B1%Trajectory(i)%cbf == B1%Trajectory(j)%cbf) then Coupled(i, j) = 1 Coupled(j, i) = 1 - write (fmiout, *) "Coupled(i,j) set to 1 because trajs have the same CBF" - write (fmiout, *) "CBFi, CBFj", B1%Trajectory(i)%cbf, B1%Trajectory(j)%cbf + !write (fmiout, *) "Coupled(i,j) set to 1 because trajs have the same CBF" + !write (fmiout, *) "CBFi, CBFj", B1%Trajectory(i)%cbf, B1%Trajectory(j)%cbf cycle end if end do @@ -717,16 +811,16 @@ subroutine FMS_BuildCoupled(B1, Coupled, selectionTime) else - do i = 1, ntraj - write (fmiout, *) "Building coupled for a triplet?", B1%Trajectory(i)%triplet - end do + !do i = 1, ntraj + ! write (fmiout, *) "Building coupled for a triplet?", B1%Trajectory(i)%triplet + !end do call FMS_BuildCoupled_ESS(B1, Coupled) - write (fmiOut, *) "Coupled matrix is:", Coupled, ", Size is:", size(Coupled) + !write (fmiOut, *) "Coupled matrix is:", Coupled, ", Size is:", size(Coupled) end if - write (fmiout, *) "---------------------------------------------------------------------------" + !write (fmiout, *) "---------------------------------------------------------------------------" end subroutine FMS_BuildCoupled @@ -740,18 +834,18 @@ subroutine FMS_BuildCoupled_ESS(B1, Coupled) ntraj = B1%NumTraj - write (fmiout, *) "Calculating coupling for trajs in bundle" - write (fmiout, *) " ----------------------------------------------------" - write (fmiout, *) "Are traj indices still fine here?" - write (fmiout, *) "Number trajs:", B1%NumTraj - write (fmiout, *) "Number CBFs:", B1%NCBFs - write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "Calculating coupling for trajs in bundle" + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "Are traj indices still fine here?" + !write (fmiout, *) "Number trajs:", B1%NumTraj + !write (fmiout, *) "Number CBFs:", B1%NCBFs + !write (fmiout, *) " ----------------------------------------------------" do i = 2, ntraj do j = 1, i - write (fmiout, *) "ID, CBF, is triplet?" - write (fmiout, *) "i", B1%Trajectory(i)%TrajID, B1%Trajectory(i)%CBF, B1%Trajectory(i)%triplet - write (fmiout, *) "j", B1%Trajectory(j)%TrajID, B1%Trajectory(j)%CBF, B1%Trajectory(j)%triplet + !write (fmiout, *) "ID, CBF, is triplet?" + !write (fmiout, *) "i", B1%Trajectory(i)%TrajID, B1%Trajectory(i)%CBF, B1%Trajectory(i)%triplet + !write (fmiout, *) "j", B1%Trajectory(j)%TrajID, B1%Trajectory(j)%CBF, B1%Trajectory(j)%triplet if (B1%Trajectory(i)%cbf == B1%Trajectory(j)%cbf) cycle if (abs(FMS_bH(B1, i, j)) > gldStochaThresh) then Coupled(i, j) = 1 @@ -885,8 +979,8 @@ subroutine FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, & integer(kind=DefInt) :: ntraj, i, j, l, itraj - write (fmiout, *) "=========================================" - write (fmiout, *) "What happens in GroupIntoBlocks?" + !write (fmiout, *) "=========================================" + !write (fmiout, *) "What happens in GroupIntoBlocks?" ntraj = B1%NumTraj @@ -901,31 +995,31 @@ subroutine FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, & if (ztrajdone(j)) cycle blocktrajid(1, i) = j ztrajdone(j) = .true. - write (fmiout, *) "The hitherto unsorted traj ", j, "was added to the block as the first traj" + !write (fmiout, *) "The hitherto unsorted traj ", j, "was added to the block as the first traj" exit end do itraj = blocktrajid(1, i) !first trajectory in block - write (fmiout, *) "This is blocktrajid inside the i loop, i = ", i - do j = 1, B1%NumTraj + 1 - write (fmiout, *) blocktrajid(j,:) - end do + !write (fmiout, *) "This is blocktrajid inside the i loop, i = ", i + !do j = 1, B1%NumTraj + 1 + ! write (fmiout, *) blocktrajid(j,:) + !end do ! If no trajectories added, then exit loop if (itraj == 0) then nblock = i - 1 - write (fmiout, *) "nblock is now ", nblock + !write (fmiout, *) "nblock is now ", nblock exit end if ! Add other coupled trajectories to this current block l = 1 do j = 1, ntraj - write (fmiout, *) "Now sorting traj j = ", j + !write (fmiout, *) "Now sorting traj j = ", j if (ztrajdone(j)) cycle !if this trajectory has already been sorted skip if (Coupled(j, itraj) == 1) then - write (fmiout, *) "traj ", j, "is in the same block as traj ", itraj + !write (fmiout, *) "traj ", j, "is in the same block as traj ", itraj l = l + 1 blocktrajid(l, i) = j ztrajdone(j) = .true. @@ -935,15 +1029,15 @@ subroutine FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, & end do - write (fmiout, *) "This is the finished blocktrajid: " - do j = 1, B1%NumTraj + 1 - write (fmiout, *) blocktrajid(j,:) - end do + !write (fmiout, *) "This is the finished blocktrajid: " + !do j = 1, B1%NumTraj + 1 + ! write (fmiout, *) blocktrajid(j,:) + !end do - write (fmiout, *) "These trajectories are in the first block: ", blocktrajid(:,1) - write (fmiout, *) "and these are in the second block: ", blocktrajid(:,2) + !write (fmiout, *) "These trajectories are in the first block: ", blocktrajid(:,1) + !write (fmiout, *) "and these are in the second block: ", blocktrajid(:,2) - write (fmiout, *) "=========================================" + !write (fmiout, *) "=========================================" end subroutine FMS_GroupIntoBlocks @@ -1019,7 +1113,7 @@ subroutine FMS_SelectBlock(ntraj, nblock, ntrajblock, totBlockpop, & ! Use random number to select which block to collapse to drand = fms_ranb(i4zero) -! find which block was randomly selected +! find which block was randomly selected ---------------- eq (14) SSAIMS paper Blockpopnormsum = 0.0d0 do i = 1, nblock Blockpopnorm = Blockpop(i) / totBlockpop @@ -1062,10 +1156,17 @@ subroutine FMS_RenormalizeBlockAmplitudes(B1, slctBlockpop, slctNtrajblock, & ! Renormalize selected block to total initial Norm (need not be 1) NormInit = FMS_Norm(B1) + !write (fmiout, *) "Inside FMS_RenormalizeBlockAmplitude: " + !write (fmiout, *) "This is NormInit", NormInit + !write (fmiout, *) "This is slctBlockpop", slctBlockpop + !write (fmiout, *) "This is slctNtrajblock", slctNtrajblock cnorm = dcmplx(sqrt(NormInit / slctBlockpop)) + !write (fmiout, *) "Renormalising the amplitudes (Inside perform_stochastic_selection): " do j = 1, slctNtrajblock jtraj = slctBlocktrajid(j) + !write (fmiout, *) "Before (jtraj, Amp): ", jtraj, B1%Trajectory(jtraj)%Amplitude B1%Trajectory(jtraj)%Amplitude = B1%Trajectory(jtraj)%Amplitude * cnorm + !write (fmiout, *) "After (jtraj, Amp): ", jtraj, B1%Trajectory(jtraj)%Amplitude end do end subroutine FMS_RenormalizeBlockAmplitudes @@ -1085,30 +1186,30 @@ subroutine FMS_KillOtherBlocks(B1, nblock, ntrajblock, iblockslct, blocktrajid) integer(kind=DefInt) :: i, j, l, jtraj ! Mark all other trajectories for death l = 0 - write (fmiout, *) "-----------------------------------------------" - write (fmiout, *) "Marking trajectories for death" + !write (fmiout, *) "-----------------------------------------------" + !write (fmiout, *) "Marking trajectories for death" do i = 1, nblock - write (fmiout, *) "Current block: ", i - write (fmiout, *) "was it selected? ", i == iblockslct + !write (fmiout, *) "Current block: ", i + !write (fmiout, *) "was it selected? ", i == iblockslct if (i == iblockslct) cycle do j = 1, ntrajblock(i) - write (fmiout, *) "We have a j! ", j + !write (fmiout, *) "We have a j! ", j l = l + 1 jtraj = blocktrajid(j, i) gliForceKill(l) = B1%Trajectory(jtraj)%TrajID - write (fmiout, *) "Trajectory will be killed (TrajID, CBF): ", B1%Trajectory(jtraj)%TrajID, B1%Trajectory(jtraj)%CBF + !write (fmiout, *) "Trajectory will be killed (TrajID, CBF): ", B1%Trajectory(jtraj)%TrajID, B1%Trajectory(jtraj)%CBF B1%Trajectory(jtraj)%Amplitude = dcmplx(0.0d0, 0.0d0) end do end do - write (fmiout, *) "Array of TrajIDs for dead/ non dead Trajs: ", gliForceKill, "size is: ", size(gliForceKill) - write (fmiout, *) "-----------------------------------------------" + !write (fmiout, *) "Array of TrajIDs for dead/ non dead Trajs: ", gliForceKill, "size is: ", size(gliForceKill) + !write (fmiout, *) "-----------------------------------------------" ! Remove any dead trajectory amplitudes do jtraj = 1, B1%NumDeadTraj B1%DeadTraj(jtraj)%Amplitude = dcmplx(0.0d0, 0.0d0) - write (fmiout, *) "Killing this trajectory: ", jtraj - write (fmiout, *) "Out of this many total dead trajs: ", B1%NumDeadTraj - write (fmiout, *) "Was it a triplet?", B1%DeadTraj(jtraj)%triplet + !write (fmiout, *) "Killing this trajectory: ", jtraj + !write (fmiout, *) "Out of this many total dead trajs: ", B1%NumDeadTraj + !write (fmiout, *) "Was it a triplet?", B1%DeadTraj(jtraj)%triplet end do end subroutine FMS_KillOtherBlocks @@ -1195,13 +1296,17 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) logical :: is_addCen ! Number of particles integer(kind=DefInt) :: npart + + ! Renormalization of the separated Bundles: + real(kind=DefReal), dimension(B1%NumStates) :: OldNorms + real(kind=DefReal) :: DNorm0 npart = B1%Trajectory(1)%NumParticles - write (fmiout, *) " ----------------------------------------------------" - write (fmiout, *) "Total number of trajs and CBFs in unseparated Bundle: " - write (fmiout, *) B1%NumTraj, B1%NCBFs - write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "Total number of trajs and CBFs in unseparated Bundle: " + !write (fmiout, *) B1%NumTraj, B1%NCBFs + !write (fmiout, *) " ----------------------------------------------------" !write (fmiout, *) " ----------------------------------------------------" !write (fmiout, *) "Number of trajs and CBFs that will be in unseparated Bundle: " @@ -1210,16 +1315,16 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) !write (fmiout, *) " ----------------------------------------------------" ! Determine number of trajs and CBFs in each StoSel block - write (fmiout, *) " ----------------------------------------------------" - write (fmiout, *) "What we currently have in the B1 bundle: " - do iTraj = 1, B1%NumTraj + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "What we currently have in the B1 bundle: " + !do iTraj = 1, B1%NumTraj - write (fmiout, *) "Info for trajectory with ID ", B1%Trajectory(iTraj)%TrajID - write (fmiout, *) "Triplet?", B1%Trajectory(iTraj)%triplet - write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms - write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF + ! write (fmiout, *) "Info for trajectory with ID ", B1%Trajectory(iTraj)%TrajID + ! write (fmiout, *) "Triplet?", B1%Trajectory(iTraj)%triplet + ! write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms + ! write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF - end do + !end do NumTrajBlock = 0 NumCBFBlock = 0 @@ -1253,15 +1358,15 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) end do - write (fmiout, *) "Number of Trajs and CBFs in block ", iblocks - write (fmiout, *) NumTrajBlock(iblocks), NumCBFBlock(iblocks) + !write (fmiout, *) "Number of Trajs and CBFs in block ", iblocks + !write (fmiout, *) NumTrajBlock(iblocks), NumCBFBlock(iblocks) end do - write (fmiout, *) "This is BlockCBFID :" - do iblocks = 1, BundlesMultDim - write (fmiout, *) BlockCBFID(:, iblocks) - end do + !write (fmiout, *) "This is BlockCBFID :" + !do iblocks = 1, BundlesMultDim + ! write (fmiout, *) BlockCBFID(:, iblocks) + !end do !Create temporary bundles, separated according 2 selection mode: ! S only, or T only, (as for S/T combined, no SS should happen @@ -1287,29 +1392,30 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) ! Copy over to the temp, separated bundles: ! - matching trajs according 2 BlockTrajID column ! - matching Centroids (make sure CentID matches) - write (fmiout, *) "The blocktrajid we are using for copying trajs to sep Bundles: " - do iTraj = 1, B1%NumTraj + 1 - write (fmiout, *) blocktrajid(iTraj,:) - end do - write (fmiout, *) " ----------------------------------------------------" - write (fmiout, *) "Original centroids in B1" - write (fmiout, *) "Does B1 have smth like positions for Centroids?", B1%Centroids(1)%Particle(1)%get_pos() - do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) - write (fmiout, *) "Centroid number ", iTraj - write (fmiout, *) "is centroid to trajectories ", B1%Centroids(iTraj)%CentID - write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() - end do + !write (fmiout, *) "The blocktrajid we are using for copying trajs to sep Bundles: " + !do iTraj = 1, B1%NumTraj + 1 + ! write (fmiout, *) blocktrajid(iTraj,:) + !end do + + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "Original centroids in B1" + !write (fmiout, *) "Does B1 have smth like positions for Centroids?", B1%Centroids(1)%Particle(1)%get_pos() + !do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) + ! write (fmiout, *) "Centroid number ", iTraj + ! write (fmiout, *) "is centroid to trajectories ", B1%Centroids(iTraj)%CentID + ! write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() + !end do - write (fmiout, *) " ----------------------------------------------------" - write (fmiout, *) "Filling up the empty, separated bundles with Trajs from B1:" + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "Filling up the empty, separated bundles with Trajs from B1:" do iblocks = 1, BundlesMultDim - write (fmiout, *) "=====================================================" - write (fmiout, *) "Now adding trajs to sep Bundle ", iblocks - write (fmiout, *) "Should start with traj of index ", blocktrajid(1,iblocks) - write (fmiout, *) "=====================================================" + !write (fmiout, *) "=====================================================" + !write (fmiout, *) "Now adding trajs to sep Bundle ", iblocks + !write (fmiout, *) "Should start with traj of index ", blocktrajid(1,iblocks) + !write (fmiout, *) "=====================================================" CenCount = 1 CBFcount = 1 @@ -1320,10 +1426,10 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) ! iTraj --> ID traj will have in temp bundle ! addTraj --> ID of corresponding traj in B1 - write (fmiout, *) "iTraj ", iTraj - write (fmiout, *) "represents traj", addTraj - write (fmiout, *) "The BlockCBFID for current block: ", BlockCBFID(:, iblocks) - write (fmiout, *) "has actual CBF", B1%Trajectory(addTraj)%CBF, "the BlockCBFID CBF is:", BlockCBFID(CBFcount,iblocks) + !write (fmiout, *) "iTraj ", iTraj + !write (fmiout, *) "represents traj", addTraj + !write (fmiout, *) "The BlockCBFID for current block: ", BlockCBFID(:, iblocks) + !write (fmiout, *) "has actual CBF", B1%Trajectory(addTraj)%CBF, "the BlockCBFID CBF is:", BlockCBFID(CBFcount,iblocks) !if (iblocks>1) then ! write (fmiout, *) "in sep Bundle, needs to have CBF", B1%Trajectory(addTraj)%CBF - NumCBFBlock(iblocks-1) @@ -1341,7 +1447,7 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) ! (BlockCBFID(CBFcount,iblocks) is the CBF of previous the iteration) if (B1%Trajectory(addTraj)%CBF /= BlockCBFID(CBFcount,iblocks)) then !CBFcurr = B1%Trajectory(iTraj)%CBF - write (fmiout, *) "Adding one to CBFcount" + !write (fmiout, *) "Adding one to CBFcount" CBFcount = CBFcount + 1 end if @@ -1358,9 +1464,9 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) do jblocks = 1, iblocks-1 BlockCBFCount = BlockCBFCount + NumCBFBlock(jblocks) end do - write (fmiout, *) "BlockCBFcount (all CBFs added to previous blocks: ", BlockCBFCount + !write (fmiout, *) "BlockCBFcount (all CBFs added to previous blocks: ", BlockCBFCount BundlesMult(iblocks)%Trajectory(iTraj)%CBF = B1%Trajectory(addTraj)%CBF - BlockCBFCount - write (fmiout, *) "CBF receives the number (must start at 1)", B1%Trajectory(addTraj)%CBF - BlockCBFCount + !write (fmiout, *) "CBF receives the number (must start at 1)", B1%Trajectory(addTraj)%CBF - BlockCBFCount else BundlesMult(iblocks)%Trajectory(iTraj)%CBF = B1%Trajectory(addTraj)%CBF end if @@ -1370,13 +1476,13 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) ! Careful here, figuring this out has been a mess but we (we=Vera) pray that it works now ! - do we really cover all possible CBF ID pairs, i.e. catch all centroids? ! --------------------------------------------------------------------------------------- - write (fmiout, *) "___________________________" + !write (fmiout, *) "___________________________" ! - start from the 1st CBF ID (i.e. 1st row, current column of BlockCBFID) ! and go over all CBF IDs as CBFcount increases CBF_i = BlockCBFID(CBFcount, iblocks) - write (fmiout, *) "CBF_i is, CBF_j will be", CBF_i, CBF_i + 1 - write (fmiout, *) "loop limit will be:", B1%NCBFs + !write (fmiout, *) "CBF_i is, CBF_j will be", CBF_i, CBF_i + 1 + !write (fmiout, *) "loop limit will be:", B1%NCBFs ! - loop over possible CBF ID pairs and check 4 matching CentID pair in B1 is_addCen = .false. @@ -1386,22 +1492,22 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) ! if max of cens in ! temp bundle is reached - write (fmiout, *) "CBF_i and CBF_j", CBF_i, CBF_j + !write (fmiout, *) "CBF_i and CBF_j", CBF_i, CBF_j if (any(CBF_i == BlockCBFID(:,iblocks)) .and. any(CBF_j == BlockCBFID(:,iblocks)) .and. CBF_i /= CBF_j) then - write (fmiout, *) "dream come true", CBF_j, CBF_i - write (fmiout, *) "CenID is", CenID + !write (fmiout, *) "dream come true", CBF_j, CBF_i + !write (fmiout, *) "CenID is", CenID ! For this CentID pair, find out corresponding ID of Centroid in B1 (CenID) do Cen_TrajID = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) - write(fmiout, *) "Is Cen ", Cen_TrajID, "the correct one?" - write (fmiout, *) B1%Centroids(Cen_TrajID)%CentID + !write(fmiout, *) "Is Cen ", Cen_TrajID, "the correct one?" + !write (fmiout, *) B1%Centroids(Cen_TrajID)%CentID Cen_i = B1%Centroids(Cen_TrajID)%CentID(1); Cen_j = B1%Centroids(Cen_TrajID)%CentID(2) if (CBF_i == Cen_i .and. CBF_j == Cen_j .or. CBF_j == Cen_i .and. CBF_i == Cen_j) then CenID = Cen_TrajID - write (fmiout, *) "The ID of the centroid we need is ", CenID + !write (fmiout, *) "The ID of the centroid we need is ", CenID is_addCen = .true. !cycle end if @@ -1413,7 +1519,7 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) BundlesMult(iblocks)%Centroids(CenCount) = B1%Centroids(CenID) BundlesMult(iblocks)%Centroids(CenCount)%CentID = [CBF_j, CBF_i] BundlesMult(iblocks)%Centroids(CenCount)%TrajID = CenCount - write (fmiout, *) "We added Cen", CenID, "of B1 to block Bundle ", iblocks, "as Cen number", CenCount + !write (fmiout, *) "We added Cen", CenID, "of B1 to block Bundle ", iblocks, "as Cen number", CenCount if (iblocks>1) then !BlockCBFCount = 0 @@ -1425,57 +1531,60 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) BundlesMult(iblocks)%Centroids(CenCount)%CentID = [CBF_j - BlockCBFCount, & CBF_i - BlockCBFCount] - write (fmiout, *) "Cen receives the number (must start at 1)" + !write (fmiout, *) "Cen receives the number (must start at 1)" end if CenCount = CenCount + 1 cycle end if - else - write (fmiout, *) "Index pair is not the one needed for the CentIDs in this block" + !else + ! write (fmiout, *) "Index pair is not the one needed for the CentIDs in this block" end if end do - write (fmiout, *) "___________________________" + !write (fmiout, *) "___________________________" end do end do - write (fmiout, *) "Do we have Centroids after copying them?", glzCentroids - do iblocks = 1, BundlesMultDim - if (BundlesMult(iblocks)%NCBFs > 1) then - do iTraj = 1, (((BundlesMult(iblocks)%NCBFs - 1) * BundlesMult(iblocks)%NCBFs) / 2) - write (fmiout, *) "Centroid number ", iTraj - write (fmiout, *) "is centroid to trajectories ", BundlesMult(iblocks)%Centroids(iTraj)%CentID - write (fmiout, *) "And has position: ", BundlesMult(iblocks)%Centroids(iTraj)%Particle(1)%get_pos() - end do - else - write (fmiout, *) "This block only has one CBF: ", iblocks, ". Therefore no centroids." - end if - end do + !write (fmiout, *) "Do we have Centroids after copying them?", glzCentroids + !do iblocks = 1, BundlesMultDim + ! if (BundlesMult(iblocks)%NCBFs > 1) then + ! do iTraj = 1, (((BundlesMult(iblocks)%NCBFs - 1) * BundlesMult(iblocks)%NCBFs) / 2) + ! write (fmiout, *) "Centroid number ", iTraj + ! write (fmiout, *) "is centroid to trajectories ", BundlesMult(iblocks)%Centroids(iTraj)%CentID + ! write (fmiout, *) "And has position: ", BundlesMult(iblocks)%Centroids(iTraj)%Particle(1)%get_pos() + ! end do + ! else + ! write (fmiout, *) "This block only has one CBF: ", iblocks, ". Therefore no centroids." + ! end if + !end do - write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) " ----------------------------------------------------" - write (fmiout, *) "What we filled into the separated bundle: " - do iblocks = 1, BundlesMultDim - do iTraj = 1, BundlesMult(iblocks)%NumTraj + !write (fmiout, *) "What we filled into the separated bundle: " + !do iblocks = 1, BundlesMultDim + ! do iTraj = 1, BundlesMult(iblocks)%NumTraj - write (fmiout, *) "Info for trajectory with ID ", BundlesMult(iblocks)%Trajectory(iTraj)%TrajID - write (fmiout, *) "Triplet?", BundlesMult(iblocks)%Trajectory(iTraj)%triplet - write (fmiout, *) "Ms", BundlesMult(iblocks)%Trajectory(iTraj)%Ms - write (fmiout, *) "CBF", BundlesMult(iblocks)%Trajectory(iTraj)%CBF + ! write (fmiout, *) "Info for trajectory with ID ", BundlesMult(iblocks)%Trajectory(iTraj)%TrajID + ! write (fmiout, *) "Triplet?", BundlesMult(iblocks)%Trajectory(iTraj)%triplet + ! write (fmiout, *) "Ms", BundlesMult(iblocks)%Trajectory(iTraj)%Ms + ! write (fmiout, *) "CBF", BundlesMult(iblocks)%Trajectory(iTraj)%CBF - end do - end do - write (fmiout, *) " ----------------------------------------------------" + ! end do + !end do + !write (fmiout, *) " ----------------------------------------------------" ! Normalise the separate bundles before SS (and keep the og norm, amplitudes?) + !DNorm0 = 1.d0 !do iblocks = 1, BundlesMultDim - ! call FMS_Renormalize(BundlesMult(i), + ! call FMS_Branching(BundlesMult(iblocks), OldNorms) ! no idea yet what OldNorms should actually + ! ! be but just to try it out + ! call FMS_Renormalize(BundlesMult(iblocks), OldNorms, DNorm0) !end do end subroutine getMultBundles @@ -1501,11 +1610,11 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, ! All we need to do is to adjust gliForceKill to contain the TrajIDs for B1 ! not those referring to the temp separated bundles - do i = 1, B1%NumTraj - write (fmiout, *) "These trajs are in the original B1 (before copying back)" - write (fmiout, *) "TrajID: ", B1%Trajectory(i)%TrajID, "DeadTime: ", B1%Trajectory(i)%DeadTime - write (fmiout, *) "IsDead? (checking if DeadTime < CurrentTime)", B1%Trajectory(i)%DeadTime < B1%CurrentTime - end do + !do i = 1, B1%NumTraj + ! write (fmiout, *) "These trajs are in the original B1 (before copying back)" + ! write (fmiout, *) "TrajID: ", B1%Trajectory(i)%TrajID, "DeadTime: ", B1%Trajectory(i)%DeadTime + ! write (fmiout, *) "IsDead? (checking if DeadTime < CurrentTime)", B1%Trajectory(i)%DeadTime < B1%CurrentTime + !end do ! How to get the information on which Trajs are dead? !write (fmiout, *) "These trajs are dead in the original B1 (before copying back)" @@ -1520,55 +1629,55 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, do iTraj = 1, B1%NumTraj ! Does this loop always account for all dead and alive trajs??? ! Even after several kills in the Bundle? ! Do we need to add B1%NumDeadTraj somehow asp? - write (fmiout, *) " ------------------------" - write (fmiout, *) "Current iTraj: ", iTraj + !write (fmiout, *) " ------------------------" + !write (fmiout, *) "Current iTraj: ", iTraj - if (any(B1%Trajectory(:)%TrajID == iTraj)) then - write(fmiout, *) "This traj is alive" - end if + !if (any(B1%Trajectory(:)%TrajID == iTraj)) then + ! write(fmiout, *) "This traj is alive" + !end if if (any(B1%DeadTraj(:)%TrajID == iTraj)) then - write(fmiout, *) "This traj is dead" + !write(fmiout, *) "This traj is dead" ! This TrajID needs to be remembered and 'left blank' when temp separated Bundles are copied back DeadCount = DeadCount + 1 DeadTrajIDs(DeadCount) = iTraj end if if (any(B1%DeadTraj(:)%TrajID == iTraj) .and. any(B1%Trajectory(:)%TrajID == iTraj)) then - write(fmiout, *) "This traj is dead and alive" + !write(fmiout, *) "This traj is dead and alive" call FMS_DieError("This shouldn't happen") end if - write (fmiout, *) " ------------------------" + !write (fmiout, *) " ------------------------" end do - write (fmiout, *) "These are our dead trajectories: ", DeadTrajIDs + !write (fmiout, *) "These are our dead trajectories: ", DeadTrajIDs hitDead = .false. do i = 1, BundlesMultDim if (StoSelMode(i) == 1 .or. StoSelMode(i) == 3) then - write (fmiout, *) "Killing happened in block Bundle ", i - write (fmiout, *) "Trajs need to be copied back to B1 to keep track of who died" - - do iTraj = 1, BundlesMult(i)%NumTraj - write (fmiout, *) " ------------------------" - write (fmiout, *) "Copying back Traj ", iTraj - write (fmiout, *) "Ms", BundlesMult(i)%Trajectory(iTraj)%Ms - write (fmiout, *) "CBF", BundlesMult(i)%Trajectory(iTraj)%CBF - write (fmiout, *) "What is the amplitude?", BundlesMult(i)%Trajectory(iTraj)%Amplitude - write (fmiout, *) " ------------------------" - write (fmiout, *) "Is this traj still the same in B1?", B1%Trajectory(iTraj)%TrajID - write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms - write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF - write (fmiout, *) "What is the amplitude?", B1%Trajectory(iTraj)%Amplitude - write (fmiout, *) " ------------------------" - end do - write (fmiout, *) " ----------------------------------------------------" - write (fmiout, *) "Now, actually copy and check again: " - write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "Killing happened in block Bundle ", i + !write (fmiout, *) "Trajs need to be copied back to B1 to keep track of who died" + + !do iTraj = 1, BundlesMult(i)%NumTraj + ! write (fmiout, *) " ------------------------" + ! write (fmiout, *) "Copying back Traj ", iTraj + ! write (fmiout, *) "Ms", BundlesMult(i)%Trajectory(iTraj)%Ms + ! write (fmiout, *) "CBF", BundlesMult(i)%Trajectory(iTraj)%CBF + ! write (fmiout, *) "What is the amplitude?", BundlesMult(i)%Trajectory(iTraj)%Amplitude + ! write (fmiout, *) " ------------------------" + ! write (fmiout, *) "Is this traj still the same in B1?", B1%Trajectory(iTraj)%TrajID + ! write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms + ! write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF + ! write (fmiout, *) "What is the amplitude?", B1%Trajectory(iTraj)%Amplitude + ! write (fmiout, *) " ------------------------" + !end do + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "Now, actually copy and check again: " + !write (fmiout, *) " ----------------------------------------------------" ! Copy back trajectories, adjust temp bundle TrajID, CBF back to refer to B1 - write (fmiout, *) "Copying trajs for block bundle ", i + !write (fmiout, *) "Copying trajs for block bundle ", i TrajCount = 0 CBFCount = 0 @@ -1577,12 +1686,13 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, ! does not represent the total number of ! blocks anymore ! can we store the Coupled(B1) and use here? + ! TODO: 27/04/26: this problem is still relevant!! do j = 1, i - 1 TrajCount = TrajCount + BundlesMult(j)%NumTraj CBFCount = CBFCount + BundlesMult(j)%NCBFs end do - write (fmiout, *) "The", i-1, "previous bundles had: " - write (fmiout, *) TrajCount, "Trajs and", CBFcount, "CBFs" + !write (fmiout, *) "The", i-1, "previous bundles had: " + !write (fmiout, *) TrajCount, "Trajs and", CBFcount, "CBFs" end if @@ -1591,48 +1701,61 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, B1%Trajectory(iTraj+TrajCount) = BundlesMult(i)%Trajectory(iTraj) !TODO: figure out what original TrajID, CBF ID was if (any(BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj == DeadTrajIDs(:))) then - write (fmiout, *) "We hit the dead traj (iTraj)", & - BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj - write (fmiout, *) "Adding the DeadCount ", DeadCount, "from now on" + !write (fmiout, *) "We hit the dead traj (iTraj)", & + ! BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + !write (fmiout, *) "Adding the DeadCount ", DeadCount, "from now on" hitDead = .true. end if if (hitDead) then B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + & DeadCount + B1%Trajectory(iTraj+TrajCount)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude ! used to break stuff else B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + B1%Trajectory(iTraj+TrajCount)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude ! used to break stuff end if B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF + BundlesMult(i-1)%NCBFs write (fmiout, *) "Traj of index ", iTraj, "was copied over and received ID, ", & B1%Trajectory(iTraj+TrajCount)%TrajID + + ! copy over amplitudes for first block as well: + ! as it's the first block we don't need to worry about: + ! - how many trajs we had in previous blocks (for sure) + ! - how many dead trajs (not too sure, but we pray for it each night) + else if (i==1) then + B1%Trajectory(iTraj)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude end if - end do - do iTraj = 1, BundlesMult(i)%NumTraj - write (fmiout, *) " ------------------------" - write (fmiout, *) "Copying back Traj ", iTraj - write (fmiout, *) "Ms", BundlesMult(i)%Trajectory(iTraj)%Ms - write (fmiout, *) "CBF", BundlesMult(i)%Trajectory(iTraj)%CBF - write (fmiout, *) "What is the amplitude?", BundlesMult(i)%Trajectory(iTraj)%Amplitude - write (fmiout, *) " ------------------------" - write (fmiout, *) "Is this traj still the same in B1?", B1%Trajectory(iTraj)%TrajID - write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms - write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF - write (fmiout, *) "What is the amplitude?", B1%Trajectory(iTraj)%Amplitude - write (fmiout, *) " ------------------------" end do + + ! 23/04/26: Definitely go back here and check that copying back is done correctly + ! also apart from amplitudes + + !do iTraj = 1, BundlesMult(i)%NumTraj + ! write (fmiout, *) " ------------------------" + ! write (fmiout, *) "Copying back Traj ", iTraj + ! write (fmiout, *) "Ms", BundlesMult(i)%Trajectory(iTraj)%Ms + ! write (fmiout, *) "CBF", BundlesMult(i)%Trajectory(iTraj)%CBF + ! write (fmiout, *) "What is the amplitude?", BundlesMult(i)%Trajectory(iTraj)%Amplitude + ! write (fmiout, *) " ------------------------" + ! write (fmiout, *) "Is this traj still the same in B1?", B1%Trajectory(iTraj)%TrajID + ! write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms + ! write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF + ! write (fmiout, *) "What is the amplitude?", B1%Trajectory(iTraj)%Amplitude + ! write (fmiout, *) " ------------------------" + !end do end if end do - write (fmiout, *) " ----------------------------------------------------" - write (fmiout, *) "These trajs are in B1 (after copying back)" - do i = 1, B1%NumTraj - write (fmiout, *) "TrajID: ", B1%Trajectory(i)%TrajID, "CBF: ", B1%Trajectory(i)%CBF - write (fmiout, *) "Ms: ", B1%Trajectory(i)%Ms - write (fmiout, *) "Amp: ", B1%Trajectory(i)%Amplitude - end do - write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) " ----------------------------------------------------" + !write (fmiout, *) "These trajs are in B1 (after copying back)" + !do i = 1, B1%NumTraj + ! write (fmiout, *) "TrajID: ", B1%Trajectory(i)%TrajID, "CBF: ", B1%Trajectory(i)%CBF + ! write (fmiout, *) "Ms: ", B1%Trajectory(i)%Ms + ! write (fmiout, *) "Amp: ", B1%Trajectory(i)%Amplitude + !end do + !write (fmiout, *) " ----------------------------------------------------" end subroutine copy_MultBundles_to_original_bundle @@ -1707,4 +1830,66 @@ subroutine copy_state_bundles_to_original_bundle(B1, BundleSS) end subroutine copy_state_bundles_to_original_bundle +! Copied over from PropagationModule for now +! TODO: figure out how to properly get this from the PropagationModule +! (by making the subroutine public in there and then adjusting +! the dependencies in the MakeFile??) + +!> +!! Renormalizes a trajectory bundle such that Branching(Bundle)=OldNorms. +!! +!! This is necessary because frozen Gaussian propagation does not +!! conserve wavefunction normalization. +!! \param OldNorms Norms of each trajectory to normalize to +!! \param DNorm0 Overall normalization at time t=0; this is usually +!! 1 except when pruning trajectories on restart. +!! @ingroup propagation +!< + subroutine FMS_Renormalize(B1, OldNorms, DNorm0) + use BundleCalcsModule, only: FMS_Norm, FMS_Branching + type(T_TrajectoryBundle), intent(inout) :: B1 + real(kind=DefReal), intent(in) :: OldNorms(B1%NumStates) + real(kind=DefReal), intent(in) :: DNorm0 + + real(kind=DefReal), allocatable, save :: NewNorms(:) + integer(kind=DefInt), save :: LastDimension + real(kind=DefReal) :: DNorm + integer(kind=DefInt) :: IState, ITraj + + !write (fmiout, *) "Inside FMS_Renormalize: " + !write (fmiout, *) "This is OldNorms: ", OldNorms + + if (B1%NumStates /= LastDimension) then + if (LastDimension /= 0) deallocate (NewNorms) + allocate (NewNorms(B1%NumStates)) + LastDimension = B1%NumStates + end if + + call FMS_Branching(B1, NewNorms) + + !write (fmiout, *) "Inside FMS_Renormalize: " + !write (fmiout, *) "This is NewNorms: ", NewNorms + + do IState = 1, B1%NumStates + if (abs(NewNorms(IState)) > FPZero) then + NewNorms(IState) = sqrt(OldNorms(IState) / NewNorms(IState)) + else + NewNorms(IState) = 0.d0 + end if + end do + + do ITraj = 1, B1%NumTraj + B1%Trajectory(ITraj)%Amplitude = NewNorms(B1%Trajectory(ITraj)%StateID) * B1%Trajectory(ITraj)%Amplitude + end do + +! Now make sure overall normalization is unity... + DNorm = FMS_Norm(B1) + DNorm = sqrt(DNorm0 / DNorm) + + do ITraj = 1, B1%NumTraj + B1%Trajectory(ITraj)%Amplitude = DNorm * B1%Trajectory(ITraj)%Amplitude + end do + + end subroutine FMS_Renormalize + end module SelectionModule diff --git a/src/remove_dead.f90 b/src/remove_dead.f90 index 70574ec..841f70d 100644 --- a/src/remove_dead.f90 +++ b/src/remove_dead.f90 @@ -97,7 +97,7 @@ subroutine FMS_RemoveDead(B1) NotCoupled = all(Coupled(:, i) == 0) OnIgnoreState = (StateID == glIgnoreState) PopBelowThresh = (Population(i) < spdPopToSpawn) -! GAIMS changed +! GAIMS changed --- using total pop of all sublevels pop = 0.d0 do n = 1, ntraj if (B1%Trajectory(n)%CBF == B1%Trajectory(i)%CBF) then @@ -110,7 +110,7 @@ subroutine FMS_RemoveDead(B1) .or. (PopBelowThresh .and. NotCoupled)) ForceDead = any(gliForceKill(:) == TrajID) if (ForceDead) then - B1%Trajectory(i)%DeadTime = -9999.0d0 + B1%Trajectory(i)%DeadTime = -9999.0d0 ! Why is this done? end if MarkForDeath = (MarkForDeath .or. ForceDead) @@ -220,6 +220,26 @@ subroutine FMS_RemoveDead(B1) B1 = BTemp call BTemp%destroy() + write (fmiout, *) "Is Bundle still okay after remove_dead reconstruction?" + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "These trajs are in B1" + do iTraj = 1, B1%NumTraj + write (fmiout, *) "TrajID: ", B1%Trajectory(iTraj)%TrajID, "CBF: ", B1%Trajectory(iTraj)%CBF + write (fmiout, *) "Ms: ", B1%Trajectory(iTraj)%Ms + write (fmiout, *) "Amp: ", B1%Trajectory(iTraj)%Amplitude + end do + write (fmiout, *) " ----------------------------------------------------" + + write (fmiout, *) "Centroids of B1" + write (fmiout, *) "Part 3: after remove_dead" + do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) + write (fmiout, *) "Centroid number", iTraj + write (fmiout, *) "Is centroid to trajectories", B1%Centroids(iTraj)%CentID + write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() + end do + write (fmiout, *) " ----------------------------------------------------" + + end if end subroutine FMS_RemoveDead From e7dcbeb133a9e62db81b126f8397e0cbd7cf335c Mon Sep 17 00:00:00 2001 From: verab98 Date: Mon, 11 May 2026 09:37:35 +0100 Subject: [PATCH 08/10] Getting closer to a working version Relevant change: - loop for translating indices which trajectory to kill now starts at 1 instead of 2 --- src/modules/SelectionModule.f90 | 94 +++++++++++++++++++++------------ src/remove_dead.f90 | 4 ++ 2 files changed, 63 insertions(+), 35 deletions(-) diff --git a/src/modules/SelectionModule.f90 b/src/modules/SelectionModule.f90 index b94c01d..b41d643 100644 --- a/src/modules/SelectionModule.f90 +++ b/src/modules/SelectionModule.f90 @@ -137,6 +137,7 @@ subroutine FMS_StochasticCollapse(B1) else if (NTrip /= 0) then NumDeadTrajSave = B1%NumDeadTraj + gliForceKill = 0 !write (fmiout, *) "Asking this question once and for all: do we call FMS_StochasticCollapse at each time step??" ! yes we do... @@ -398,16 +399,16 @@ subroutine FMS_StochasticCollapse(B1) SaveForceKill = gliForceKill - do iTraj = 2, B1%NumTraj + 1 + do iTraj = 1, B1%NumTraj + 1 ! Why did I ever start this one from 2?? do j = 1, B1%NumTraj + 1 if (MergedSepBlockTrajID(iTraj,j) == 0) cycle do DeadID = 1, size(SaveForceKill) !write (fmiout, *) "many numbers (?)", SaveForceKill(DeadID) if (MergedSepBlockTrajID(iTraj,j) == SaveForceKill(DeadID)) then - !write (fmiout, *) "Match for gliForceKill: ", SaveForceKill(DeadID), "(this should be 4)" - !write (fmiout, *) "with MergedSepblock: ", MergedSepBlockTrajID(iTraj,j) - !write (fmiout, *) "i and j are: ", iTraj, j, "(should be 2, 4)" - !write (fmiout, *) "Corresponding MergedBlock: ", MergedBlockTrajID(iTraj,j), "(this should be 6)" + write (fmiout, *) "Match for gliForceKill: ", SaveForceKill(DeadID) + write (fmiout, *) "with MergedSepblock: ", MergedSepBlockTrajID(iTraj,j) + write (fmiout, *) "i and j are: ", iTraj, j + write (fmiout, *) "Corresponding MergedBlock: ", MergedBlockTrajID(iTraj,j) gliForceKill(DeadID) = MergedBlockTrajID(iTraj,j) end if end do @@ -418,6 +419,11 @@ subroutine FMS_StochasticCollapse(B1) write (fmiout, *) gliForceKill end if + if (IsSelection) then + write (fmiout, *) "(Before copying back to original Bundle)" + write (fmiout, *) "Number of Dead Trajs is now: ", B1%NumDeadTraj + end if + !if (B1%NumDeadTraj > 0) then ! do iTraj = 1, B1%NumDeadTraj @@ -494,28 +500,35 @@ subroutine FMS_StochasticCollapse(B1) ! StoSel (but not sure if that breaks things somewhere...?) - call copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode) - !write (fmiout, *) "States were copied back without error" + if (IsSelection) then + call copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode) + write (fmiout, *) "States were copied back without error" + end if - write (fmiout, *) "After copy_MultBundles_to_original_Bundle, the B1 norm is: " - write (fmiout, *) "Norm: ", FMS_Norm(B1) - write (fmiout, *) "and the corresponding amplitudes are" - do iTraj = 1, B1%NumTraj - write (fmiout, *) "Trajectory ", iTraj, "AMplitude: ", B1%Trajectory(iTraj)%Amplitude - end do - write (fmiout, *) "-----------------------------------------------------------" + if (IsSelection) then + write (fmiout, *) "(After copying back to original Bundle)" + write (fmiout, *) "Number of Dead Trajs is now: ", B1%NumDeadTraj - !end if + write (fmiout, *) "After copy_MultBundles_to_original_Bundle, the B1 norm is: " + write (fmiout, *) "Norm: ", FMS_Norm(B1) + write (fmiout, *) "and the corresponding amplitudes are" + do iTraj = 1, B1%NumTraj + write (fmiout, *) "Trajectory ", iTraj, "has ID", B1%Trajectory(iTraj)%TrajID, "and AMplitude: ", & + B1%Trajectory(iTraj)%Amplitude + end do + end if + + write (fmiout, *) "-----------------------------------------------------------" - !write (fmiout, *) " ----------------------------------------------------" - !write (fmiout, *) "Check that Centroids of B1 remain unaffected" - !write (fmiout, *) "Part 2: B1 Centroids after" - !do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) - ! write (fmiout, *) "Centroid number", iTraj - ! write (fmiout, *) "Is centroid to trajectories", B1%Centroids(iTraj)%CentID - ! write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() - !end do - !write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "Check that Centroids of B1 remain unaffected" + write (fmiout, *) "Part 2: B1 Centroids after" + do iTraj = 1, (((B1%NCBFs - 1) * B1%NCBFs) / 2) + write (fmiout, *) "Centroid number", iTraj + write (fmiout, *) "Is centroid to trajectories", B1%Centroids(iTraj)%CentID + write (fmiout, *) "And has position: ", B1%Centroids(iTraj)%Particle(1)%get_pos() + end do + write (fmiout, *) " ----------------------------------------------------" ! Destroy the temporary block bundles do i = 1, nblockMerged @@ -523,6 +536,14 @@ subroutine FMS_StochasticCollapse(B1) end do deallocate(SaveNorms, SaveTrajAmps, BundlesMult) + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "Check that TrajIDs of B1 remain unaffected" + write (fmiout, *) "(At the end of StochasticCollapse)" + do iTraj = 1, B1%NumTraj + write (fmiout, *) "Traj number", iTraj, "has TrajID", B1%Trajectory(iTraj)%TrajID + end do + write (fmiout, *) " ----------------------------------------------------" + !call FMS_DieError("Now, please go back and reconsider your life choices!") !else @@ -1701,19 +1722,22 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, B1%Trajectory(iTraj+TrajCount) = BundlesMult(i)%Trajectory(iTraj) !TODO: figure out what original TrajID, CBF ID was if (any(BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj == DeadTrajIDs(:))) then - !write (fmiout, *) "We hit the dead traj (iTraj)", & - ! BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj - !write (fmiout, *) "Adding the DeadCount ", DeadCount, "from now on" + write (fmiout, *) "We hit the dead traj (iTraj)", & + BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + write (fmiout, *) "Adding the DeadCount ", DeadCount, "from now on" hitDead = .true. end if if (hitDead) then B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + & DeadCount + !B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF + ! TODO: 30/04/26: Do we have to adjust CBF IDs just like for TrajID?? B1%Trajectory(iTraj+TrajCount)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude ! used to break stuff else B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj B1%Trajectory(iTraj+TrajCount)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude ! used to break stuff + !B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF end if B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF + BundlesMult(i-1)%NCBFs write (fmiout, *) "Traj of index ", iTraj, "was copied over and received ID, ", & @@ -1748,14 +1772,14 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, end if end do - !write (fmiout, *) " ----------------------------------------------------" - !write (fmiout, *) "These trajs are in B1 (after copying back)" - !do i = 1, B1%NumTraj - ! write (fmiout, *) "TrajID: ", B1%Trajectory(i)%TrajID, "CBF: ", B1%Trajectory(i)%CBF - ! write (fmiout, *) "Ms: ", B1%Trajectory(i)%Ms - ! write (fmiout, *) "Amp: ", B1%Trajectory(i)%Amplitude - !end do - !write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "These trajs are in B1 (after copying back)" + do i = 1, B1%NumTraj + write (fmiout, *) "TrajID: ", B1%Trajectory(i)%TrajID, "CBF: ", B1%Trajectory(i)%CBF + write (fmiout, *) "Ms: ", B1%Trajectory(i)%Ms + write (fmiout, *) "Amp: ", B1%Trajectory(i)%Amplitude + end do + write (fmiout, *) " ----------------------------------------------------" end subroutine copy_MultBundles_to_original_bundle diff --git a/src/remove_dead.f90 b/src/remove_dead.f90 index 841f70d..503fa6f 100644 --- a/src/remove_dead.f90 +++ b/src/remove_dead.f90 @@ -218,6 +218,10 @@ subroutine FMS_RemoveDead(B1) BTemp%NCBFs = nCBF B1 = BTemp + + write (fmiout, *) "(After remove_dead)" + write (fmiout, *) "Number of Dead Trajs is now: ", B1%NumDeadTraj + call BTemp%destroy() write (fmiout, *) "Is Bundle still okay after remove_dead reconstruction?" From 9e3f374a73ec1a917c3eb6e7374bb237fdb2c6d4 Mon Sep 17 00:00:00 2001 From: verab98 Date: Wed, 13 May 2026 10:17:19 +0100 Subject: [PATCH 09/10] More fixes fixed the reordering of blocktrajid into blocks that correspond to the stochastic selection modes (i.e. all Singlet-only block trajids in one column, all Triplet-only in one column, before this did not work if a block was spread over more than two individual blocks) Also, activated a set of useful prints for debugging --- src/modules/SelectionModule.f90 | 97 +++++++++++++++++++++++++++------ 1 file changed, 80 insertions(+), 17 deletions(-) diff --git a/src/modules/SelectionModule.f90 b/src/modules/SelectionModule.f90 index b41d643..c66a545 100644 --- a/src/modules/SelectionModule.f90 +++ b/src/modules/SelectionModule.f90 @@ -80,7 +80,8 @@ subroutine FMS_StochasticCollapse(B1) logical :: isSing, isTrip, isSingTrip logical :: IsSelection ! --- Added 4 sep SS end --- - integer(kind=DefInt) :: iTraj, i, j, k, nblock_outer, totDeadCount + integer(kind=DefInt) :: iTraj, i, j, k, nblock_outer, totDeadCount, prev_add_traj + integer(kind=DefInt) :: skip_positions ! AIMSWISS: Shall we perform stochastic selection logical :: performSelection ! AIMSWISS: Current selection time @@ -149,10 +150,10 @@ subroutine FMS_StochasticCollapse(B1) !write (fmiout, *) "Building coupled matrix for total Bundle (B1, unseparated)" Coupled = 0 call FMS_BuildCoupled(B1, Coupled, selectionTime) - !write (fmiout, *) "This is the total, converged coupled matrix of size ", size(Coupled) - !do iTraj = 1, B1%NumTraj - ! write (fmiout, *) Coupled(iTraj,:) - !end do + write (fmiout, *) "This is the total, converged coupled matrix of size ", size(Coupled) + do iTraj = 1, B1%NumTraj + write (fmiout, *) Coupled(iTraj,:) + end do !write (fmiout, *) "Get number of blocks for total Bundle (B1, unseparated)" blocktrajid = 0 @@ -212,19 +213,59 @@ subroutine FMS_StochasticCollapse(B1) ! Merge blocks depending on isSIng, isTrip ! (No action for isSingTrip because not to be touched in StoSel + write (fmiout, *) "blocktrajid before merging: " + do i = 1, B1%NumTraj + 1 + write (fmiout, *) blocktrajid(i,:) + end do + do i = 1, nblock_outer + !write (fmiout, *) "Entering loop for i = ", i + !write (fmiout, *) "Part of blocktrajid we are using: ", blocktrajid(i,:) + !write (fmiout, *) "General question: what is NumTrajBlock???", NumTrajBlock + if (StoSelMode(i) == 0) cycle + !prev_add_traj = 0 + do j = i+1, nblock_outer + write (fmiout, *) "Entering loop for j = ", j if (StoSelMode(i) == StoSelMode(j)) then + ! Copy over block j TrajIDs to i column of blocktrajid - blocktrajid(NumTrajBlock(i)+1:NumTrajBlock(i)+NumTrajBlock(j), i) = blocktrajid(1:NumTrajBlock(j),j) + ! Set the NumTrajBlock(i:j) positions in the i column to the indices in the jth column + !blocktrajid(NumTrajBlock(i)+1:NumTrajBlock(i)+NumTrajBlock(j), i) = blocktrajid(1:NumTrajBlock(j),j) + !write (fmiout, *) "Copied over ", blocktrajid(1:NumTrajBlock(j),j), " to ",& + ! blocktrajid( sum(NumTrajBlock(1:j)) : NumTrajBlock(i) + NumTrajBlock(j), i) + + ! Count how many positions need to be skipped + skip_positions = 0 + do k = 1, j-1 + if (StoSelMode(k) == StoSelMode(j)) then + skip_positions = skip_positions + NumTrajBlock(k) + end if + end do + !write (fmiout, *) "This many positions ", skip_positions, "need 2 be skipped in col", i + + prev_add_traj = 0 + !write (fmiout, *) "++++++++++++++++++++" + do k = 1, NumTrajBlock(j) + !write (fmiout, *) "Number of trajs added to column", i, "in a previous step: ", prev_add_traj + !write (fmiout, *) "Number of trajs in previous blocks: ", NumTrajBlock(1:j-1), "in sum: ", & + ! sum(NumTrajBlock(1:j-1)) + !write (fmiout, *) "Now replace element ", skip_positions+k, "in ", & + ! i, "th column by", blocktrajid(k,j) + blocktrajid(skip_positions+k, i) = blocktrajid(k,j) + prev_add_traj = prev_add_traj + 1 + end do + !write (fmiout, *) "++++++++++++++++++++" + ! Set j column to zero blocktrajid(:,j) = 0 - ! Set StoSelMode to zero for block copied over (the j column) - StoSelMode(j) = 0 + + !write (fmiout, *) "Part of blocktrajid was modified to: ", blocktrajid(i,:) + !write (fmiout, *) "The ith column is now: ", blocktrajid(:,i) end if @@ -232,6 +273,19 @@ subroutine FMS_StochasticCollapse(B1) end do + ! Set StoSelMode to zero if column of blocktrajid is fully zero + ! i.e. if the block was copied over (the j column) + do i = 1, nblock_outer + if (all(blocktrajid(:,i) == 0)) then + StoSelMode(i) = 0 + end if + end do + + write (fmiout, *) "blocktrajid after merging: " + do i = 1, B1%NumTraj + 1 + write (fmiout, *) blocktrajid(i,:) + end do + MergedBlockTrajID = blocktrajid nblockMerged = nblock_outer ! TODO: For now, nblockMerged is not used, but do we need to account for change ! in number of blocks after merging (so, nblockMerged < nblock_outer depending @@ -314,6 +368,7 @@ subroutine FMS_StochasticCollapse(B1) end if end do nblockMerged = BundlesMultDim + write (fmiout, *) "nblock_outer = ", nblock_outer, "StoSelMode = ", StoSelMode !write (fmiout, *) "BundlesMultDim based on merged blocktrajid: ", BundlesMultDim !write (fmiout, *) "nblock based on merged blocktrajid: ", nblockMerged @@ -323,6 +378,7 @@ subroutine FMS_StochasticCollapse(B1) ! Create BundlesMultDim many temp bundles and fill them up ! with corresponding trajs from B1 + write (fmiout, *) "This is the blocktrajid to be used in getMultBundles: ", blocktrajid call getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) ! Set the norm of each temporary bundle to one (for SS) @@ -1364,8 +1420,13 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) do iblocks = 1, bundlesmultdim CBFcount = 0 + !CBFcurr = 1 + CBFCount do iTraj =1, B1%NumTraj + !write (fmiout, *) "Counting NCBFs in block ", iblocks + !write (fmiout, *) "The current CBF is: ", CBFcurr + !write (fmiout, *) "This is blocktrajid", blocktrajid + ! count traj if its ID is in the current block if (any(blocktrajid(:, iblocks) == iTraj)) then NumTrajBlock(iblocks) = NumTrajBlock(iblocks) + 1 @@ -1379,8 +1440,9 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) end do - !write (fmiout, *) "Number of Trajs and CBFs in block ", iblocks - !write (fmiout, *) NumTrajBlock(iblocks), NumCBFBlock(iblocks) + ! useful printing --> keep! + write (fmiout, *) "Number of Trajs and CBFs in block ", iblocks + write (fmiout, *) NumTrajBlock(iblocks), NumCBFBlock(iblocks) end do @@ -1402,13 +1464,14 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) end do - !write (fmiout, *) " ----------------------------------------------------" - !write (fmiout, *) "These are the empty, separated Bundles we created: " - !do iblocks = 1, BundlesMultDim - ! write (fmiout, *) "Block number ", iblocks - ! write (fmiout, *) "Number trajs, Number CBFs:", BundlesMult(iblocks)%NumTraj, BundlesMult(iblocks)%NCBFs - !end do - !write (fmiout, *) " ----------------------------------------------------" + ! useful printing --> keep! + write (fmiout, *) " ----------------------------------------------------" + write (fmiout, *) "These are the empty, separated Bundles we created: " + do iblocks = 1, BundlesMultDim + write (fmiout, *) "Block number ", iblocks + write (fmiout, *) "Number trajs, Number CBFs:", BundlesMult(iblocks)%NumTraj, BundlesMult(iblocks)%NCBFs + end do + write (fmiout, *) " ----------------------------------------------------" ! Copy over to the temp, separated bundles: ! - matching trajs according 2 BlockTrajID column From e00ba9019d57ad4b5bd231293bd976ba053930dd Mon Sep 17 00:00:00 2001 From: verab98 Date: Wed, 20 May 2026 18:01:38 +0100 Subject: [PATCH 10/10] Is this a working version? - realised that the TrajIDs can remain unchanged when copying from B1 into the spin-separated bundles - realised that, with the TrajIDs remaining unchanged, all that needs to be copied back to B1 after a selection are the amplitudes --- src/modules/SamplingModule.f90 | 4 +- src/modules/SelectionModule.f90 | 460 +++++++++++++++++++------------- src/openfms.F90 | 2 +- 3 files changed, 285 insertions(+), 181 deletions(-) diff --git a/src/modules/SamplingModule.f90 b/src/modules/SamplingModule.f90 index 5946f68..9ea1d6d 100644 --- a/src/modules/SamplingModule.f90 +++ b/src/modules/SamplingModule.f90 @@ -1371,7 +1371,9 @@ subroutine FMS_InitialSwarm(B1) B1%Trajectory(n)%TrajID = n call B1%Trajectory(n)%set_pos(X + dX) call B1%Trajectory(n)%set_mom(P + dP) ! useless here because dP is zero, but for completeness... - dX = [10.0d0, 0.d0, 0.d0] + dX = [10.0d0, 0.d0, 0.d0] ! used for placing two Gaussians in front of, one after + ! the crossing point + !dX = [-0.25d0, 0.d0, 0.d0] ! for having three Gaussians in front of crossing end do else call B1%Trajectory(ntraj)%set_pos(X + dX) diff --git a/src/modules/SelectionModule.f90 b/src/modules/SelectionModule.f90 index c66a545..95bfc26 100644 --- a/src/modules/SelectionModule.f90 +++ b/src/modules/SelectionModule.f90 @@ -405,7 +405,8 @@ subroutine FMS_StochasticCollapse(B1) write (fmiout, *) "Norm: ", FMS_Norm(B1) write (fmiout, *) "and the corresponding amplitudes are" do iTraj = 1, B1%NumTraj - write (fmiout, *) "Trajectory ", iTraj, "AMplitude: ", B1%Trajectory(iTraj)%Amplitude + write (fmiout, *) "Trajectory ", iTraj, "TrajID ", B1%Trajectory(iTraj)%TrajID, & + "AMplitude: ", B1%Trajectory(iTraj)%Amplitude end do write (fmiout, *) "----------------------------" @@ -429,52 +430,62 @@ subroutine FMS_StochasticCollapse(B1) if (StoSelMode(i) == 1 .or. StoSelMode(i) == 3) then write (fmiout, *) "StoSelMode for block ", i, "is ", StoSelMode(i) + IsSelection = .false. ! Do the stochastic selection if (BundlesMult(i)%NCBFs > 1 ) then write (fmiout, *) "Performing stochastic selection for block Bundle", i, & "with", BundlesMult(i)%NCBFs, "CBFs" - IsSelection = .false. call perform_stochastic_selection(BundlesMult(i), selectionTime, IsSelection) ! Adjust gliForceKill only if selection actually happened - if (i>1 .and. IsSelection) then - !write (fmiout, *) "Adjusting gliForceKill using MergedSepBlockTrajID (this is the full matrix): " - !do j = 1, B1%NumTraj + 1 - ! write (fmiout, *) MergedSepBlockTrajID(j,:) - !end do - - !write (fmiout, *) "and MergedBlockTrajID (this is the full matrix): " - !do j = 1, B1%NumTraj + 1 - ! write (fmiout, *) MergedBlockTrajID(j,:) - !end do - - write (fmiout, *) "gliForcKill for block Bundle", i - write (fmiout, *) gliForceKill - !write (fmiout, *) "adjusting to B1 indiced..." - - SaveForceKill = gliForceKill - - do iTraj = 1, B1%NumTraj + 1 ! Why did I ever start this one from 2?? - do j = 1, B1%NumTraj + 1 - if (MergedSepBlockTrajID(iTraj,j) == 0) cycle - do DeadID = 1, size(SaveForceKill) - !write (fmiout, *) "many numbers (?)", SaveForceKill(DeadID) - if (MergedSepBlockTrajID(iTraj,j) == SaveForceKill(DeadID)) then - write (fmiout, *) "Match for gliForceKill: ", SaveForceKill(DeadID) - write (fmiout, *) "with MergedSepblock: ", MergedSepBlockTrajID(iTraj,j) - write (fmiout, *) "i and j are: ", iTraj, j - write (fmiout, *) "Corresponding MergedBlock: ", MergedBlockTrajID(iTraj,j) - gliForceKill(DeadID) = MergedBlockTrajID(iTraj,j) - end if - end do - end do + if (IsSelection) then + + write (fmiout, *) "Number of dead trajs: ", B1%NumDeadTraj + do j = 1, B1%NumDeadTraj + write (fmiout, *) "TrajID of dead traj", j, ": ", B1%DeadTraj(j)%TrajID end do - write (fmiout, *) "Adjusted gliForcKill for block Bundle", i, "and obtained" - write (fmiout, *) gliForceKill end if + ! Adjust gliForceKill only if selection actually happened + !if (i>1 .and. IsSelection) then + ! !write (fmiout, *) "Adjusting gliForceKill using MergedSepBlockTrajID (this is the full matrix): " + ! !do j = 1, B1%NumTraj + 1 + ! ! write (fmiout, *) MergedSepBlockTrajID(j,:) + ! !end do + + ! !write (fmiout, *) "and MergedBlockTrajID (this is the full matrix): " + ! !do j = 1, B1%NumTraj + 1 + ! ! write (fmiout, *) MergedBlockTrajID(j,:) + ! !end do + + ! write (fmiout, *) "gliForcKill for block Bundle", i + ! write (fmiout, *) gliForceKill + ! !write (fmiout, *) "adjusting to B1 indiced..." + + ! SaveForceKill = gliForceKill + + ! do iTraj = 1, B1%NumTraj + 1 ! Why did I ever start this one from 2?? + ! do j = 1, B1%NumTraj + 1 + ! if (MergedSepBlockTrajID(iTraj,j) == 0) cycle + ! do DeadID = 1, size(SaveForceKill) + ! !write (fmiout, *) "many numbers (?)", SaveForceKill(DeadID) + ! if (MergedSepBlockTrajID(iTraj,j) == SaveForceKill(DeadID)) then + ! write (fmiout, *) "Match for gliForceKill: ", SaveForceKill(DeadID) + ! write (fmiout, *) "with MergedSepblock: ", MergedSepBlockTrajID(iTraj,j) + ! write (fmiout, *) "i and j are: ", iTraj, j + ! write (fmiout, *) "Corresponding MergedBlock: ", MergedBlockTrajID(iTraj,j) + ! gliForceKill(DeadID) = MergedBlockTrajID(iTraj,j) + ! end if + ! end do + ! end do + ! end do + + ! write (fmiout, *) "Adjusted gliForcKill for block Bundle", i, "and obtained" + ! write (fmiout, *) gliForceKill + !end if + if (IsSelection) then write (fmiout, *) "(Before copying back to original Bundle)" write (fmiout, *) "Number of Dead Trajs is now: ", B1%NumDeadTraj @@ -539,7 +550,8 @@ subroutine FMS_StochasticCollapse(B1) do i = 1, BundlesMultDim write (fmiout, *) "Block: ", i do iTraj = 1, BundlesMult(i)%NumTraj - write (fmiout, *) "Trajectory ", iTraj, "AMplitude: ", BundlesMult(i)%Trajectory(iTraj)%Amplitude + write (fmiout, *) "Trajectory ", iTraj, "TrajID ", BundlesMult(i)%Trajectory(iTraj)%TrajID, & + "AMplitude: ", BundlesMult(i)%Trajectory(iTraj)%Amplitude end do end do write (fmiout, *) "----------------------------" @@ -556,8 +568,25 @@ subroutine FMS_StochasticCollapse(B1) ! StoSel (but not sure if that breaks things somewhere...?) + write (fmiout, *) "blocktrajid right before copy_MultBundles_to_B1: " + do i = 1, B1%NumTraj + 1 + write (fmiout, *) blocktrajid(i,:) + end do + + write (fmiout, *) "MergedBlockTrajID: " + do i = 1, B1%NumTraj + 1 + write (fmiout, *) MergedBlockTrajID(i,:) + end do + + write (fmiout, *) "MergedSepBlockTrajID: " + do i = 1, B1%NumTraj + 1 + write (fmiout, *) MergedBlockTrajID(i,:) + end do + if (IsSelection) then - call copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode) + call copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode, & + MergedBlockTrajID, MergedSepBlockTrajID) + ! 18/05/26: Don't we need to normalise here ???!!?!?! write (fmiout, *) "States were copied back without error" end if @@ -664,7 +693,9 @@ subroutine perform_stochastic_selection(B1, selectionTime, IsSelection) ntrajblock = 0 nblock = 0 call FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, ntrajblock, nblock) + write (fmiout, *) "Right before performing SS:" write (fmiout, *) "Number of Blocks: ", nblock + write (fmiout, *) "Coupled is: ", Coupled ! There should be at least one block of trajectories if (nblock == 0) then @@ -1106,10 +1137,10 @@ subroutine FMS_GroupIntoBlocks(B1, Coupled, blocktrajid, & end do - !write (fmiout, *) "This is the finished blocktrajid: " - !do j = 1, B1%NumTraj + 1 - ! write (fmiout, *) blocktrajid(j,:) - !end do + write (fmiout, *) "This is the finished blocktrajid: " + do j = 1, B1%NumTraj + 1 + write (fmiout, *) blocktrajid(j,:) + end do !write (fmiout, *) "These trajectories are in the first block: ", blocktrajid(:,1) !write (fmiout, *) "and these are in the second block: ", blocktrajid(:,2) @@ -1263,8 +1294,8 @@ subroutine FMS_KillOtherBlocks(B1, nblock, ntrajblock, iblockslct, blocktrajid) integer(kind=DefInt) :: i, j, l, jtraj ! Mark all other trajectories for death l = 0 - !write (fmiout, *) "-----------------------------------------------" - !write (fmiout, *) "Marking trajectories for death" + write (fmiout, *) "-----------------------------------------------" + write (fmiout, *) "Marking trajectories for death" do i = 1, nblock !write (fmiout, *) "Current block: ", i !write (fmiout, *) "was it selected? ", i == iblockslct @@ -1274,12 +1305,12 @@ subroutine FMS_KillOtherBlocks(B1, nblock, ntrajblock, iblockslct, blocktrajid) l = l + 1 jtraj = blocktrajid(j, i) gliForceKill(l) = B1%Trajectory(jtraj)%TrajID - !write (fmiout, *) "Trajectory will be killed (TrajID, CBF): ", B1%Trajectory(jtraj)%TrajID, B1%Trajectory(jtraj)%CBF + write (fmiout, *) "Trajectory will be killed (TrajID, CBF): ", B1%Trajectory(jtraj)%TrajID, B1%Trajectory(jtraj)%CBF B1%Trajectory(jtraj)%Amplitude = dcmplx(0.0d0, 0.0d0) end do end do - !write (fmiout, *) "Array of TrajIDs for dead/ non dead Trajs: ", gliForceKill, "size is: ", size(gliForceKill) - !write (fmiout, *) "-----------------------------------------------" + write (fmiout, *) "Array of TrajIDs for dead/ non dead Trajs: ", gliForceKill, "size is: ", size(gliForceKill) + write (fmiout, *) "-----------------------------------------------" ! Remove any dead trajectory amplitudes do jtraj = 1, B1%NumDeadTraj @@ -1510,10 +1541,10 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) ! iTraj --> ID traj will have in temp bundle ! addTraj --> ID of corresponding traj in B1 - !write (fmiout, *) "iTraj ", iTraj - !write (fmiout, *) "represents traj", addTraj - !write (fmiout, *) "The BlockCBFID for current block: ", BlockCBFID(:, iblocks) - !write (fmiout, *) "has actual CBF", B1%Trajectory(addTraj)%CBF, "the BlockCBFID CBF is:", BlockCBFID(CBFcount,iblocks) + write (fmiout, *) "iTraj ", iTraj + write (fmiout, *) "represents traj", addTraj + write (fmiout, *) "The BlockCBFID for current block: ", BlockCBFID(:, iblocks) + write (fmiout, *) "has actual CBF", B1%Trajectory(addTraj)%CBF, "the BlockCBFID CBF is:", BlockCBFID(CBFcount,iblocks) !if (iblocks>1) then ! write (fmiout, *) "in sep Bundle, needs to have CBF", B1%Trajectory(addTraj)%CBF - NumCBFBlock(iblocks-1) @@ -1539,7 +1570,9 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) ! set TrajID to iTraj to make sure indices in temp bundle start at 1 call BundlesMult(iblocks)%Trajectory(iTraj)%create(npart, B1%NumStates) BundlesMult(iblocks)%Trajectory(iTraj) = B1%Trajectory(addTraj) - BundlesMult(iblocks)%Trajectory(iTraj)%TrajID = iTraj + !BundlesMult(iblocks)%Trajectory(iTraj)%TrajID = iTraj + !BundlesMult(iblocks)%Trajectory(iTraj)%TrajID = addTraj + BundlesMult(iblocks)%Trajectory(iTraj)%TrajID = B1%Trajectory(addTraj)%TrajID ! check if CBF IDs need to be reset to 1 (the case if we are not in the first separated block) ! to make sure CBF IDs also start at 1 @@ -1649,18 +1682,18 @@ subroutine getMultBundles(B1, BundlesMult, BundlesMultDim, blocktrajid) !write (fmiout, *) " ----------------------------------------------------" - !write (fmiout, *) "What we filled into the separated bundle: " - !do iblocks = 1, BundlesMultDim - ! do iTraj = 1, BundlesMult(iblocks)%NumTraj + write (fmiout, *) "What we filled into the separated bundle: " + do iblocks = 1, BundlesMultDim + do iTraj = 1, BundlesMult(iblocks)%NumTraj - ! write (fmiout, *) "Info for trajectory with ID ", BundlesMult(iblocks)%Trajectory(iTraj)%TrajID - ! write (fmiout, *) "Triplet?", BundlesMult(iblocks)%Trajectory(iTraj)%triplet - ! write (fmiout, *) "Ms", BundlesMult(iblocks)%Trajectory(iTraj)%Ms - ! write (fmiout, *) "CBF", BundlesMult(iblocks)%Trajectory(iTraj)%CBF + write (fmiout, *) "Info for trajectory with ID ", BundlesMult(iblocks)%Trajectory(iTraj)%TrajID + write (fmiout, *) "Triplet?", BundlesMult(iblocks)%Trajectory(iTraj)%triplet + write (fmiout, *) "Ms", BundlesMult(iblocks)%Trajectory(iTraj)%Ms + write (fmiout, *) "CBF", BundlesMult(iblocks)%Trajectory(iTraj)%CBF - ! end do - !end do - !write (fmiout, *) " ----------------------------------------------------" + end do + end do + write (fmiout, *) " ----------------------------------------------------" ! Normalise the separate bundles before SS (and keep the og norm, amplitudes?) @@ -1676,16 +1709,18 @@ end subroutine getMultBundles ! Copy over trajectories after the selection from BundlesMult ! to the original bundle. Note that we must copy over ! also the trajectories marked for death. - subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode) + subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, StoSelMode, & + MergedBlockTrajID, MergedSepBlockTrajID) integer(kind=DefInt), intent(in) :: BundlesMultDim integer(kind=DefInt), dimension(BundlesMultDim), intent(in) :: StoSelMode type(T_TrajectoryBundle), intent(inout) :: B1 type(T_TrajectoryBundle), intent(in) :: BundlesMult(BundlesMultDim) - !type(T_TrajectoryBundle), pointer :: BSS_i - !type(T_Trajectory), pointer :: T_i - integer(kind=DefInt) :: i, j, iTraj + integer(kind=DefInt), dimension(B1%NumTraj + 1, B1%NumTraj + 1), intent(in) :: & + MergedBlockTrajID, MergedSepBlockTrajID + integer(kind=DefInt) :: i, j, iTraj, B1TrajID, sepBTrajID integer(kind=DefInt) :: TrajCount, CBFCount, DeadCount integer(kind=DefInt), dimension(50) :: DeadTrajIDs + integer(kind=DefInt) :: add_DeadTraj, iDeadTraj logical :: hitDead !iSingTraj, iTripTraj, iCBF !integer(kind=DefInt) :: TrajIDSing, CBFIDSing, TrajIDTrip, CBFIDTrip @@ -1708,133 +1743,200 @@ subroutine copy_MultBundles_to_original_bundle(B1, BundlesMult, BundlesMultDim, ! Try a different approach to what is (or will be) commented out further down - DeadCount = 0 - DeadTrajIDs = 0 - do iTraj = 1, B1%NumTraj ! Does this loop always account for all dead and alive trajs??? - ! Even after several kills in the Bundle? - ! Do we need to add B1%NumDeadTraj somehow asp? - !write (fmiout, *) " ------------------------" - !write (fmiout, *) "Current iTraj: ", iTraj - - !if (any(B1%Trajectory(:)%TrajID == iTraj)) then - ! write(fmiout, *) "This traj is alive" - !end if - - if (any(B1%DeadTraj(:)%TrajID == iTraj)) then - !write(fmiout, *) "This traj is dead" - ! This TrajID needs to be remembered and 'left blank' when temp separated Bundles are copied back - DeadCount = DeadCount + 1 - DeadTrajIDs(DeadCount) = iTraj - end if + do i = 1, BundlesMultDim - if (any(B1%DeadTraj(:)%TrajID == iTraj) .and. any(B1%Trajectory(:)%TrajID == iTraj)) then - !write(fmiout, *) "This traj is dead and alive" - call FMS_DieError("This shouldn't happen") - end if + write (fmiout, *) "i (the how manyth block)= ", i + do j = 1, BundlesMult(i)%NumTraj - !write (fmiout, *) " ------------------------" - end do - !write (fmiout, *) "These are our dead trajectories: ", DeadTrajIDs + write (fmiout, *) "j (the how manyth traj)= ", j + write(fmiout, *) "TrajID in sepB: ", BundlesMult(i)%Trajectory(j)%TrajID - hitDead = .false. + do iTraj = 1, B1%NumTraj + B1%NumDeadTraj - do i = 1, BundlesMultDim - if (StoSelMode(i) == 1 .or. StoSelMode(i) == 3) then - !write (fmiout, *) "Killing happened in block Bundle ", i - !write (fmiout, *) "Trajs need to be copied back to B1 to keep track of who died" - - !do iTraj = 1, BundlesMult(i)%NumTraj - ! write (fmiout, *) " ------------------------" - ! write (fmiout, *) "Copying back Traj ", iTraj - ! write (fmiout, *) "Ms", BundlesMult(i)%Trajectory(iTraj)%Ms - ! write (fmiout, *) "CBF", BundlesMult(i)%Trajectory(iTraj)%CBF - ! write (fmiout, *) "What is the amplitude?", BundlesMult(i)%Trajectory(iTraj)%Amplitude - ! write (fmiout, *) " ------------------------" - ! write (fmiout, *) "Is this traj still the same in B1?", B1%Trajectory(iTraj)%TrajID - ! write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms - ! write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF - ! write (fmiout, *) "What is the amplitude?", B1%Trajectory(iTraj)%Amplitude - ! write (fmiout, *) " ------------------------" - !end do - !write (fmiout, *) " ----------------------------------------------------" - !write (fmiout, *) "Now, actually copy and check again: " - !write (fmiout, *) " ----------------------------------------------------" - - ! Copy back trajectories, adjust temp bundle TrajID, CBF back to refer to B1 - !write (fmiout, *) "Copying trajs for block bundle ", i - - TrajCount = 0 - CBFCount = 0 - if (i>1) then ! Problem here: we reduce BundlesMultDim - ! when merging Coupled, so BundlesMultDim - ! does not represent the total number of - ! blocks anymore - ! can we store the Coupled(B1) and use here? - ! TODO: 27/04/26: this problem is still relevant!! - do j = 1, i - 1 - TrajCount = TrajCount + BundlesMult(j)%NumTraj - CBFCount = CBFCount + BundlesMult(j)%NCBFs - end do - !write (fmiout, *) "The", i-1, "previous bundles had: " - !write (fmiout, *) TrajCount, "Trajs and", CBFcount, "CBFs" + add_DeadTraj = 0 - end if + if (BundlesMult(i)%Trajectory(j)%TrajID == iTraj) then + write (fmiout, *) "It's a match! for ", BundlesMult(i)%Trajectory(j)%TrajID, & + "(sepB TrajID)", "and", iTraj, "(B1 TrajID)" - do iTraj = 1, BundlesMult(i)%NumTraj - if (i>1) then - B1%Trajectory(iTraj+TrajCount) = BundlesMult(i)%Trajectory(iTraj) !TODO: figure out what original TrajID, CBF ID was - - if (any(BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj == DeadTrajIDs(:))) then - write (fmiout, *) "We hit the dead traj (iTraj)", & - BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj - write (fmiout, *) "Adding the DeadCount ", DeadCount, "from now on" - hitDead = .true. - end if + do iDeadTraj = 1, B1%NumDeadTraj + if (B1%DeadTraj(iDeadTraj)%TrajID < iTraj) then + add_DeadTraj = add_DeadTraj + 1 + end if + end do - if (hitDead) then - B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + & - DeadCount - !B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF - ! TODO: 30/04/26: Do we have to adjust CBF IDs just like for TrajID?? - B1%Trajectory(iTraj+TrajCount)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude ! used to break stuff - else - B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj - B1%Trajectory(iTraj+TrajCount)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude ! used to break stuff - !B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF - end if - B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF + BundlesMult(i-1)%NCBFs - write (fmiout, *) "Traj of index ", iTraj, "was copied over and received ID, ", & - B1%Trajectory(iTraj+TrajCount)%TrajID - - ! copy over amplitudes for first block as well: - ! as it's the first block we don't need to worry about: - ! - how many trajs we had in previous blocks (for sure) - ! - how many dead trajs (not too sure, but we pray for it each night) - else if (i==1) then - B1%Trajectory(iTraj)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude + write (fmiout, *) "The amplitude ", B1%Trajectory(iTraj - add_DeadTraj)%Amplitude, & + "will be overwritten to", BundlesMult(i)%Trajectory(j)%Amplitude + B1%Trajectory(iTraj - add_DeadTraj)%Amplitude = BundlesMult(i)%Trajectory(j)%Amplitude end if end do - ! 23/04/26: Definitely go back here and check that copying back is done correctly - ! also apart from amplitudes - - !do iTraj = 1, BundlesMult(i)%NumTraj - ! write (fmiout, *) " ------------------------" - ! write (fmiout, *) "Copying back Traj ", iTraj - ! write (fmiout, *) "Ms", BundlesMult(i)%Trajectory(iTraj)%Ms - ! write (fmiout, *) "CBF", BundlesMult(i)%Trajectory(iTraj)%CBF - ! write (fmiout, *) "What is the amplitude?", BundlesMult(i)%Trajectory(iTraj)%Amplitude - ! write (fmiout, *) " ------------------------" - ! write (fmiout, *) "Is this traj still the same in B1?", B1%Trajectory(iTraj)%TrajID - ! write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms - ! write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF - ! write (fmiout, *) "What is the amplitude?", B1%Trajectory(iTraj)%Amplitude - ! write (fmiout, *) " ------------------------" - !end do - end if + end do + end do + !DeadCount = 0 + !DeadTrajIDs = 0 + !do iTraj = 1, B1%NumTraj ! Does this loop always account for all dead and alive trajs??? + ! ! Even after several kills in the Bundle? + ! ! Do we need to add B1%NumDeadTraj somehow asp? + ! !write (fmiout, *) " ------------------------" + ! !write (fmiout, *) "Current iTraj: ", iTraj + + ! !if (any(B1%Trajectory(:)%TrajID == iTraj)) then + ! ! write(fmiout, *) "This traj is alive" + ! !end if + + ! if (any(B1%DeadTraj(:)%TrajID == iTraj)) then + ! !write(fmiout, *) "This traj is dead" + ! ! This TrajID needs to be remembered and 'left blank' when temp separated Bundles are copied back + ! DeadCount = DeadCount + 1 + ! DeadTrajIDs(DeadCount) = iTraj + ! end if + + ! if (any(B1%DeadTraj(:)%TrajID == iTraj) .and. any(B1%Trajectory(:)%TrajID == iTraj)) then + ! !write(fmiout, *) "This traj is dead and alive" + ! call FMS_DieError("Traj is dead and alive. This shouldn't happen") + ! end if + + ! !write (fmiout, *) " ------------------------" + !end do + !!write (fmiout, *) "These are our dead trajectories: ", DeadTrajIDs + + !hitDead = .false. + + !do i = 1, BundlesMultDim + ! if (StoSelMode(i) == 1 .or. StoSelMode(i) == 3) then + ! !write (fmiout, *) "Killing happened in block Bundle ", i + ! !write (fmiout, *) "Trajs need to be copied back to B1 to keep track of who died" + + ! do iTraj = 1, BundlesMult(i)%NumTraj + ! ! write (fmiout, *) " ------------------------" + ! ! write (fmiout, *) "Copying back Traj ", iTraj + ! ! write (fmiout, *) "Ms", BundlesMult(i)%Trajectory(iTraj)%Ms + ! ! write (fmiout, *) "CBF", BundlesMult(i)%Trajectory(iTraj)%CBF + ! ! write (fmiout, *) "What is the amplitude?", BundlesMult(i)%Trajectory(iTraj)%Amplitude + ! write (fmiout, *) " ------------------------" + ! write (fmiout, *) "Is this traj still the same in B1?", B1%Trajectory(iTraj)%TrajID + ! write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms + ! write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF + ! write (fmiout, *) "What is the amplitude?", B1%Trajectory(iTraj)%Amplitude + ! write (fmiout, *) " ------------------------" + ! end do + ! !write (fmiout, *) " ----------------------------------------------------" + ! !write (fmiout, *) "Now, actually copy and check again: " + ! !write (fmiout, *) " ----------------------------------------------------" + + ! ! Copy back trajectories, adjust temp bundle TrajID, CBF back to refer to B1 + ! !write (fmiout, *) "Copying trajs for block bundle ", i + ! + ! ! 20/05/26: Try something new ----------------------------------------------------- + + ! write (fmiout, *) "i (the how manyth block)= ", i + + ! do j = 1, BundlesMult(i)%NumTraj + + ! write (fmiout, *) "j (the how manyth traj)= ", j + ! add_DeadTraj = 0 + + ! do iDeadTraj = 1, DeadCount + + ! if (MergedBlockTrajID(j,i) == DeadTrajIDs(iDeadTraj)) then + ! add_DeadTraj = iDeadTraj + ! write (fmiout, *) "Hit the dead traj ", DeadTrajIDs(iDeadTraj) + ! end if + + ! end do + + ! write (fmiout, *) "add_DeadTraj is ", add_DeadTraj + + ! !B1TrajID = MergedBlockTrajID(j,i)+add_DeadTraj + ! !sepBTrajID = MergedSepBlockTrajID(j,i) + ! !write (fmiout, *) "Are these really the same as TrajIDs? -----------" + ! !write (fmiout, *) "B1TrajID", B1TrajID, "actual TrajID", & + ! ! B1%Trajectory(MergedBlockTrajID(j,i)+add_DeadTraj)%TrajID + ! !write (fmiout, *) "sepBTrajID", sepBTrajID, "actual TrajID", & + ! ! BundlesMult(i)%Trajectory(MergedSepBlockTrajID(j,i))%TrajID + ! !write (fmiout, *) "--------------" + + ! B1TrajID = B1%Trajectory(MergedBlockTrajID(j,i)+add_DeadTraj)%TrajID + ! sepBTrajID = BundlesMult(i)%Trajectory(MergedSepBlockTrajID(j,i))%TrajID + ! write (fmiout, *) "Are these really the same as TrajIDs? -----------" + ! write (fmiout, *) "B1TrajID", B1TrajID, "actual TrajID", & + ! B1%Trajectory(MergedBlockTrajID(j,i)+add_DeadTraj)%TrajID + ! write (fmiout, *) "sepBTrajID", sepBTrajID, "actual TrajID", & + ! BundlesMult(i)%Trajectory(MergedSepBlockTrajID(j,i))%TrajID + ! write (fmiout, *) "--------------" + + ! write (fmiout, *) "Attempting to copy amplitude of sep Bundle traj ", sepBTrajID, & + ! "over to traj ", B1TrajID, "in B1" + ! !write (fmiout, *) "The corresponding TrajIDs are: ", B1%Trajectory(MergedBlockTrajID(j,i)+add_DeadTraj)%TrajID, & + ! ! "in B1 and ", BundlesMult(i)%Trajectory(MergedSepBlockTrajID(j,i))%TrajID, "in sep Bundle" + ! write (fmiout, *) "The amplitude ", B1%Trajectory(B1TrajID)%Amplitude, & + ! "will be overwritten to", BundlesMult(i)%Trajectory(sepBTrajID)%Amplitude + ! B1%Trajectory(B1TrajID)%Amplitude = & + ! BundlesMult(i)%Trajectory(sepBTrajID)%Amplitude + + ! end do + + ! ! 20/05/26 end -------------------------------------------------------------------- + + ! !do iTraj = 1, BundlesMult(i)%NumTraj + ! ! if (i>1) then + ! ! B1%Trajectory(iTraj+TrajCount) = BundlesMult(i)%Trajectory(iTraj) !TODO: figure out what original TrajID, CBF ID was + + ! ! if (any(BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj == DeadTrajIDs(:))) then + ! ! write (fmiout, *) "We hit the dead traj (iTraj)", & + ! ! BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + ! ! write (fmiout, *) "Adding the DeadCount ", DeadCount, "from now on" + ! ! hitDead = .true. + ! ! end if + + ! ! if (hitDead) then + ! ! B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + & + ! ! DeadCount + ! ! !B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF + ! ! ! TODO: 30/04/26: Do we have to adjust CBF IDs just like for TrajID?? + ! ! B1%Trajectory(iTraj+TrajCount)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude ! used to break stuff + ! ! else + ! ! B1%Trajectory(iTraj+TrajCount)%TrajID = BundlesMult(i)%Trajectory(iTraj)%TrajID + BundlesMult(i-1)%NumTraj + ! ! B1%Trajectory(iTraj+TrajCount)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude ! used to break stuff + ! ! !B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF + ! ! end if + ! ! B1%Trajectory(iTraj+TrajCount)%CBF = BundlesMult(i)%Trajectory(iTraj)%CBF + BundlesMult(i-1)%NCBFs + ! ! write (fmiout, *) "Traj of index ", iTraj, "was copied over and received ID, ", & + ! ! B1%Trajectory(iTraj+TrajCount)%TrajID + + ! ! ! copy over amplitudes for first block as well: + ! ! ! as it's the first block we don't need to worry about: + ! ! ! - how many trajs we had in previous blocks (for sure) + ! ! ! - how many dead trajs (not too sure, but we pray for it each night) + ! ! else if (i==1) then + ! ! B1%Trajectory(iTraj)%Amplitude = BundlesMult(i)%Trajectory(iTraj)%Amplitude + ! ! end if + + ! !end do + + ! ! 23/04/26: Definitely go back here and check that copying back is done correctly + ! ! also apart from amplitudes + + ! do iTraj = 1, BundlesMult(i)%NumTraj + ! !write (fmiout, *) " ------------------------" + ! !write (fmiout, *) "Copying back Traj ", iTraj + ! !write (fmiout, *) "Ms", BundlesMult(i)%Trajectory(iTraj)%Ms + ! !write (fmiout, *) "CBF", BundlesMult(i)%Trajectory(iTraj)%CBF + ! !write (fmiout, *) "What is the amplitude?", BundlesMult(i)%Trajectory(iTraj)%Amplitude + ! write (fmiout, *) " ------------------------" + ! write (fmiout, *) "Is this traj still the same in B1?", B1%Trajectory(iTraj)%TrajID + ! write (fmiout, *) "Ms", B1%Trajectory(iTraj)%Ms + ! write (fmiout, *) "CBF", B1%Trajectory(iTraj)%CBF + ! write (fmiout, *) "What is the amplitude?", B1%Trajectory(iTraj)%Amplitude + ! write (fmiout, *) " ------------------------" + ! end do + ! end if + !end do + write (fmiout, *) " ----------------------------------------------------" write (fmiout, *) "These trajs are in B1 (after copying back)" do i = 1, B1%NumTraj diff --git a/src/openfms.F90 b/src/openfms.F90 index e7bb90b..4704898 100644 --- a/src/openfms.F90 +++ b/src/openfms.F90 @@ -122,7 +122,7 @@ program OpenFMS end if write (fmiOut, '(a,F0.3,a)') 'Propagate until t = ', SimulationTime, ' a.u.' - flush (fmiOut) + !flush (fmiOut) do while (Bundle%CurrentTime < SimulationTime)