Skip to content

(WiP) Separated Stochastic Selection #28

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
verab98 wants to merge 17 commits into from
Closed

(WiP) Separated Stochastic Selection #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
83a912b
Avoid initialising on a triplet state
verab98 Mar 16, 2026
af4151e
Fixing initial Amp with more than one initial TBF
verab98 Mar 16, 2026
cf5cec5
Useful print statements
verab98 Mar 16, 2026
396af33
More useful print statements
verab98 Mar 16, 2026
c2a7a00
(WiP) Separated Stochastic Selection
verab98 Apr 7, 2026
50d6fe7
Separated StoSel (still WiP)
verab98 Apr 22, 2026
caecbe7
What looks like a version that is working
verab98 Apr 29, 2026
e7dcbeb
Getting closer to a working version
verab98 May 11, 2026
9e3f374
More fixes
verab98 May 13, 2026
e00ba90
Is this a working version?
verab98 May 20, 2026
b899b3b
Cleaning up SelectionModule
verab98 May 27, 2026
9c43c67
Clean up continued
verab98 May 27, 2026
afaa8b0
Clean up continued even further
verab98 Jun 1, 2026
3993240
Cleanup even way further
verab98 Jun 2, 2026
2dc03cb
Further clean up and added test
verab98 Jun 4, 2026
e0b21d7
Merge branch 'main' into NumInitBasis
danielhollas Jun 5, 2026
a872a5e
Update Checkpoint.txt.ref
danielhollas Jun 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 24 additions & 0 deletions src/dynamics.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -80,9 +81,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
Expand Down
38 changes: 33 additions & 5 deletions src/modules/BundleCalcsModule.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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))

Expand Down
7 changes: 6 additions & 1 deletion src/modules/OverlapModule.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
! Copyright Todd J. Martinez and Raphael D. Levine, 1994
module OverlapModule
use GlobalModule
use ParticleModule
Expand Down Expand Up @@ -630,6 +629,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

Expand Down Expand Up @@ -676,6 +678,9 @@ 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, *) "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
call FMS_DieError('overlap_V_trajectory :: centroids dont match')
Expand Down
8 changes: 5 additions & 3 deletions src/modules/PropagationModule.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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))
Expand Down
97 changes: 71 additions & 26 deletions src/modules/SamplingModule.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1338,8 +1338,11 @@ subroutine FMS_InitialSwarm(B1)
real(DefReal), parameter :: SIGMA = 0.1d0 ! size of random perturbation
real(defReal) :: P_norm ! norm of momentum
real(defReal) :: B1_norm
real(defReal) :: dspl ! displacement in x direction of X, for tests w/ SOC model
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
Expand All @@ -1355,18 +1358,52 @@ 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
!dspl = 0.1
dspl = -0.25
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] ! used for placing two Gaussians in front of, one after
! the crossing point
! the setting for running the SSGAIMS test
!dX = [-0.25d0, 0.d0, 0.d0] ! for having three Gaussians in front of crossing
! ------------------------------------------------------------------------------------------
!dX = dX + [dspl, 0.d0, 0.d0] ! for having many Gaussians floating around
!dspl = dspl + 0.1 ! somewhere and somehow
!dX = dX + [dspl, 0.d0, 0.d0] ! Used for comp with full GAIMS
!dspl = dspl - 0.25
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
Expand All @@ -1379,36 +1416,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
Loading
Loading