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

Draft
verab98 wants to merge 5 commits into
base: main
Choose a base branch
from
Draft

(WiP) Separated Stochastic Selection #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 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
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 @@ -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
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 @@ -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

Expand Down Expand Up @@ -672,6 +674,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
86 changes: 60 additions & 26 deletions src/modules/SamplingModule.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Loading
Loading