Skip to content

Initial SPA1 Implementation for 1D Model#7

Open
adamjharrison wants to merge 8 commits intomainfrom
spa1
Open

Initial SPA1 Implementation for 1D Model#7
adamjharrison wants to merge 8 commits intomainfrom
spa1

Conversation

@adamjharrison
Copy link
Collaborator

Initial implementation of SPA1 interstate and intrastate terms for 1D Model (G. Granucci, M. Persico, G. Spighi, J. Chem. Phys., 2012, 137, 22A501)

@codecov-commenter
Copy link

codecov-commenter commented Jan 28, 2026

Codecov Report

❌ Patch coverage is 4.87805% with 39 lines in your changes missing coverage. Please review.
✅ Project coverage is 53.81%. Comparing base (e1be80b) to head (901a246).

Files with missing lines Patch % Lines
src/modules/OverlapModule.f90 0.00% 37 Missing and 2 partials ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main       #7      +/-   ##
==========================================
- Coverage   54.13%   53.81%   -0.32%     
==========================================
  Files          42       42              
  Lines        5980     6021      +41     
  Branches      800      809       +9     
==========================================
+ Hits         3237     3240       +3     
- Misses       2395     2431      +36     
- Partials      348      350       +2     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@danielhollas
Copy link
Member

danielhollas commented Jan 28, 2026

Huh, interesting, the AIMSWISS test is failing on the debug build with the following stack trace.
https://github.com/ispg-group/openfms/actions/runs/21448470092/job/61770647666?pr=7

At line 181 of file SelectionModule.f90
Fortran runtime error: Index '6' of dimension 1 of array 'blockpop' above upper bound of 3

Error termination. Backtrace:
#0  0x7fea74823e59 in ???
#1  0x7fea74824a71 in ???
#2  0x7fea74825082 in ???
#3  0x55fcb49103a5 in perform_stochastic_selection
	at /home/runner/work/openfms/openfms/src/modules/SelectionModule.f90:181
#4  0x55fcb4945aa9 in __selectionmodule_MOD_fms_stochasticcollapse
	at /home/runner/work/openfms/openfms/src/modules/SelectionModule.f90:111
#5  0x55fcb4a38a73 in fms_dynamics_
	at /home/runner/work/openfms/openfms/src/dynamics.f90:81
#6  0x55fcb4766c63 in openfms
	at /home/runner/work/openfms/openfms/src/openfms.F90:137
#7  0x55fcb4766d07 in main
	at /home/runner/work/openfms/openfms/src/openfms.F90:12

@adamjharrison when you get a chance, can you look at the code and understand how it could be influenced by your changes?

EDIT: Don't sweat about it too much, it's likely it's an pre-existing bug that's somehow exposed here?

@adamjharrison adamjharrison marked this pull request as ready for review February 6, 2026 14:27
Copy link
Member

@danielhollas danielhollas left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@adamjharrison thanks, this looks great!

I mostly have questions about naming, maybe that's something we can discuss in person next week.

I haven't verified all the various formulas, as I trust you did a good job there. But it would be great if you could update the PR description with a detailed explanation of how you tested everything.

It might be great actually to write some unit tests for this, we can discuss how to do it next week.

Comment on lines +826 to +838
! write(fmiOut,*) '########FFFF#############'
! write(fmiOut,*) 'd',T_i%TrajID,T_j%TrajID,T_i%StateID,T_j%StateID
! write(fmiOut,*) 'Ms', Msi,Msj
! write(fmiOut,*) 'S_j',S_ij
! write(fmiOut,*) 'x_i', T_i%Particle(1)%get_pos(1)
! write(fmiOut,*) 'x_j', T_j%Particle(1)%get_pos(1)
! write(fmiOut,*) 'p_i', T_i%Particle(1)%get_mom(1)
! write(fmiOut,*) 'p_j', T_j%Particle(1)%get_mom(1)
! write(fmiOut,*) 'a_i', T_i%Particle(1)%width
! write(fmiOut,*) 'a_j', T_j%Particle(1)%width
! write(fmiOut,*) 'phase_i', T_i%Phase
! write(fmiOut,*) 'phase_j', T_j%Phase
! write (fmiOut, *) '#########################'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Before we merge this we should get rid of all the debug printing to keep the code clean (can be done as a last step)

Suggested change
! write(fmiOut,*) '########FFFF#############'
! write(fmiOut,*) 'd',T_i%TrajID,T_j%TrajID,T_i%StateID,T_j%StateID
! write(fmiOut,*) 'Ms', Msi,Msj
! write(fmiOut,*) 'S_j',S_ij
! write(fmiOut,*) 'x_i', T_i%Particle(1)%get_pos(1)
! write(fmiOut,*) 'x_j', T_j%Particle(1)%get_pos(1)
! write(fmiOut,*) 'p_i', T_i%Particle(1)%get_mom(1)
! write(fmiOut,*) 'p_j', T_j%Particle(1)%get_mom(1)
! write(fmiOut,*) 'a_i', T_i%Particle(1)%width
! write(fmiOut,*) 'a_j', T_j%Particle(1)%width
! write(fmiOut,*) 'phase_i', T_i%Phase
! write(fmiOut,*) 'phase_j', T_j%Phase
! write (fmiOut, *) '#########################'

! Global parameter for the GAIMS model in ToyModelModule
real(kind=DefReal) :: glGrsigma

logical :: glzSPA
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should brainstorm a more descriptive name for this variable
(same goes for the SPA keyword)


end function nuc_dip_particlexf

! Calculates the Saddle Point Approximation 1st order term
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please expand this comment a bit more, include any citations (even those that are mentioned below).
Most notably, make sure to mention that this is an analytical computation that is only valid for a specific model. Also please include a formula of what is actually being computed here (e.g. what is the SPA1 term look like).

end function nuc_dip_particlexf

! Calculates the Saddle Point Approximation 1st order term
function SPA_1(T1, T2, S_ij) result(PotEn)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, we should have a more specific name for this I think

real(kind=DefReal) :: x_cent, x_1, x_2, p_1, p_2, sigma_G, alpha_1, alpha_2
complex(kind=DefComp) :: PotEn, dE, dSOC, roe

if (GlIMethod /= 4) then
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if (GlIMethod /= 4) then
if (glIMethod /= 4) then

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants