Baeck-An couplings for SH #186
Conversation
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #186 +/- ##
==========================================
+ Coverage 92.76% 92.83% +0.07%
==========================================
Files 47 47
Lines 6693 6759 +66
Branches 757 763 +6
==========================================
+ Hits 6209 6275 +66
Misses 472 472
Partials 12 12
Flags with carried forward coverage won't be shown. Click here to find out more.
|
|
@danielhollas could you please help me with the GFortran CI / Basic build (7) (pull_request) test failing all the time? I have all the tests passed locally and the other versions pass here too. I'm not sure what is the exact problem there. Otherwise it should be working Baeck-An implementation which I'm currently testing in Prague. |
|
This is the important bit
This likely means a bug in the code which is revealed by adding ABIN/.github/workflows/gfortran.yml Line 71 in 094ba89 Try adding
🎉 |
src/files.F90
Outdated
| open (UPOP, file=chfiles(UPOP), access=chaccess, action='write') | ||
| open (UPROB, file=chfiles(UPROB), access=chaccess, action='write') | ||
| open (UPES, file=chfiles(UPES), access=chaccess, action='write') | ||
| open (UNACME, file=chfiles(UNACME), access=chaccess, action='write') | ||
| !note, UNACME is now opened in mod_sh since we don't use it for inac==1 | ||
| open (UDOTPROD, file=chfiles(UDOTPROD), access=chaccess, action='write') |
There was a problem hiding this comment.
Let's put all of this (including debug files below) into a separate function sh_files_init() or something like that?
There was a problem hiding this comment.
Done in the latest commit.
| @@ -92,8 +96,8 @@ module mod_sh | |||
| ! for ethylene 2-state SA3 dynamics | |||
| integer :: ignore_state = 0 | |||
|
|
|||
| namelist /sh/ istate_init, nstate, substep, deltae, integ, inac, nohop, phase, decoh_alpha, popthr, ignore_state, & | |||
| nac_accu1, nac_accu2, popsumthr, energydifthr, energydriftthr, adjmom, revmom, & | |||
| namelist /sh/ istate_init, nstate, substep, deltae, integ, couplings, nohop, phase, decoh_alpha, popthr, ignore_state, & | |||
There was a problem hiding this comment.
Thanks for introducing the new keywords!
I think we should keep the old variables here for backwards compatibility (so that we don't break existing input files).
Can you also update sample_inputs/input.in.sh
There was a problem hiding this comment.
Can you also update sample_inputs/input.in.sh
Done!
About keeping inac and adjmom, if I keep them in the namelist, someone can set them but they will still be reassigned by couplings and mom_adjust, which have their default values. I'm afraid this would cause more problems. I've never seen input where these two were used so I wouldn't be worried about backwards compatibility of things done in Prague. When you changed the name for decoh_alpha, it was not backwards compatible and still people quickly adapted (again at least here in Prague). So I would stick just to the new keywords, yet I leave the executive decision up to you. :)
There was a problem hiding this comment.
What you could do is to set the couplings keyword to some default value, and if the user doesn't change it in the input, don't overwrite the old keywords. I think I already do that for ipimd and mdtype. But I don't feel too strongly.
There was a problem hiding this comment.
Just coming back to this, I agree it's okay to not support the old keywords and not complicate the code. We need to also document all of this in the documentation, but that's a separate conversation. :-)
There was a problem hiding this comment.
Just a note for myself: I don't know why I put the sh namelist at the module level, since it is not public, it would be better to put it into sh_init function, to better separate variables that are only needed there (such as the new couplings keyword.
Thanks for the advice. When I added it again compiled without any problems. The flag |
|
Did you run the test with the new build? These are runtime checks, not compile-time checks. You also want to disable optimizations with -O0. Only gcc7 on CI has these flags, have a look at the yaml file definition. |
Yes, I ran the tests and everything passed. I tried again also with -O0, everything passed fine. I just can't get the error locally. I will at least try to go through the last commit if there is anything problematic. |
|
What Gfortran version are you using? Can you try the following flags and |
gcc version 10.2.1 20210110 (Debian 10.2.1-6)
This worked, thanks a lot! Now it should be all working. |
|
Great! It's very weird that adding |
|
Baeck-An (BA) couplings implemented in ABIN. The implementation was tested and benchmarked on the analytic NaI model available in ABIN, see the following figure for the comparison with 2500 trajectories. BA-FSSH calculates the time-derivative coupling with the Baeck-An formula while FSSH uses the standard |
|
Awesome work! I'll have a look sometimes this week (I am currently in Switzerland on AiiDA coding week so not sure when I get to it). It would be interesting to see how Landau works on the NAI model. |
sample_inputs/input.in.sh
Outdated
| @@ -22,6 +22,8 @@ inose=0, ! NVE ensemble (thermostat turned OFF) | |||
| &sh | |||
| istate_init=2, ! initial electronic state (1 is ground state) | |||
| nstate=3, ! number of electronic states | |||
| couplings='analytic', ! non-adiabatic coupling terms 'analytic', 'baeck-an', 'none' | |||
| mom_adjust='nac' ! momentum adjustment along either 'nac' or 'velocity' | |||
There was a problem hiding this comment.
I presume this is a velocity rescaling either along NAC or by scaling the whole vector? Mom immediately associates maximum overlap for me so it might be better to pick more suitable name but I don't insist on it. (vel_rescale=nac/all?)
There was a problem hiding this comment.
Yes, it's the velocity rescaling. I see the point, we can rename it to vel_rescale, I will just wait for Dan's opinion so that I'm not changing it twice. Thanks!
There was a problem hiding this comment.
I agree that vel_rescale is better, we could even do velocity_rescaling while we're at it, it's not like we're paying for the extra characters. :-) And it's not something that users will type often (and internally in ABIN we can use a shorter name. It's good anyway to decouple the names in the input file from the names in the code, the current design is not ideally.
Perhaps a more interesting question is how to name the values.
- for the isotropic velocity rescaling, using all available kinetic energy, I would perhaps use
isotropic. I am not aware of any widely used term for this, but perhaps it would be good to see how Mario calls it (for example in Recommendations for Velocity Adjustment in Surface Hopping). - for the rescaling along NAC, the current implementation is somewhat subtle -- we first try to rescale along NAC, and if that fails, we rescale isotropically. So I think the option should be named as
nac_then_isotropicor something like that, to make it obvious.
The nac_then_isotropic approach made sense at the time (ten years ago) to reduce the number frustrated hops. It was also what Newton-X did at the time by default (I am not sure what they do now). However, as Mario recently shown quite conclusively, this can lead to artifacts for larger molecules, which will always have enough total kinetic energy to hop back. So as a follow-up to this PR, it would be great to implement the true nac option, and perhaps consider it as a safer default.
There was a problem hiding this comment.
Thanks for the comment, velocity_rescaling sounds good, I will use it. Regarding the options, I'm not a big fan of isotropic since I always saw in papers "rescaling along the momentum/velocity vector", which is also the case in Mario's paper you mentioned. Thus, I would go with velocity and nac_then_velocity, but I don't have a strong opinion.
Regarding the rescaling with NAC and then possibly velocity, this is interesting and I believe we should add also pure NAC rescaling or at least document well the current version. I'm certain most of the papers from Prague claim rescaling along NAC not knowing there is velocity rescaling in case NAC fails.
Following this topic, it might be worth making a new release of ABIN since we are all citing the old version. There has been quite a bit of development since then.
There was a problem hiding this comment.
Following this topic, it might be worth making a new release of ABIN since we are all citing the old version.
yep, it's been on my mind for some time, but was never a priority. Should we target end of the year for new release? :-) I'll create a separate issue for discussion.
There has been quite a bit of development since then.
And it's awesome to see that most of the new feature come from other people and not me.
I'm certain most of the papers from Prague claim rescaling along NAC not knowing there is velocity rescaling in case NAC fails.
Yeah, that's on me for not documenting things properly. I am pretty sure that in at least one paper I did explicitly wrote it down, but that's not an excuse.
I'm not a big fan of isotropic since I always saw in papers "rescaling along the momentum/velocity vector", which is also the case in Mario's paper you mentioned. Thus, I would go with velocity and nac_then_velocity, but I don't have a strong opinion.
Makes sense. I guess it mostly just looks a bit strange to write velocity_rescaling = 'velocity' but it does make conceptual sense in terms of "NAC vector" and "velocity vector".
Thus, I would go with velocity and nac_then_velocity, but I don't have a strong opinion.
Let's go with your suggestion and if we find something better we can change it before the next release.
There was a problem hiding this comment.
Should we target end of the year for new release? :-) I'll create a separate issue for discussion.
Sounds good, I'm up for that.
|
With the last commit and by merging #194 , all the issues raised during the review are resolved. Once you confirm @danielhollas you agree with merging, I will do so. |
|
@JanosJiri apologies for stalling this, I wanted to have a closer look at the implementation but did not yet find the time. I am a bit overwhelmed with other stuff this week so I might only get to this next week. If you want this merged go ahead and we can follow-up with a new PR if needed. Thanks for all the cool work on this! |
src/surfacehop.F90
Outdated
| @@ -165,6 +177,9 @@ subroutine sh_init(x, y, z, vx, vy, vz) | |||
| dum_fx = 0.0D0; dum_fy = 0.0D0; dum_fz = 0.0D0 | |||
| call force_clas(dum_fx, dum_fy, dum_fz, x, y, z, dum_eclas, pot) | |||
|
|
|||
| ! open nacme_all.dat for all but baeck-an couplings | |||
| if (inac /= 1) call nacmefile_init() | |||
There was a problem hiding this comment.
This should also not be opened if the couplings='none' right?
There was a problem hiding this comment.
Yes, I agree but it was opened even with couplings='none' in the original implementations so I kept it. In the latest commit, I modify it. All the tests passed locally. Hence, it shouldn't break anything in the code.
| @@ -1012,6 +1108,42 @@ subroutine interpolate(vx, vy, vz, vx_old, vy_old, vz_old, vx_int, vy_int, vz_in | |||
|
|
|||
| end subroutine interpolate | |||
|
|
|||
| ! interpolation of time-derivative coupling calculated via Baeck-An approximation | |||
| ! this routine interpolates sigma_ba between integration steps | |||
| subroutine interpolate_ba(vx, vy, vz, vx_old, vy_old, vz_old, vx_int, vy_int, vz_int, & | |||
There was a problem hiding this comment.
I am not stoked about yet another interpolation routine. At some point it would be good to refactor them into a more general routine that could then be called on various arrays. But not in this PR anyway.
| call interpolate(vx, vy, vz, vx_old, vy_old, vz_old, vx_newint, vy_newint, vz_newint, & | ||
| nacx_newint, nacy_newint, nacz_newint, en_array_newint, & | ||
| dotproduct_newint, fr) | ||
| if ((inac == 0) .or. (inac == 2)) then |
There was a problem hiding this comment.
It's not clear to me if we need to interpolate if inac==2? But better keep the existing behaviour for now.
There was a problem hiding this comment.
I don't think it is necessary to interpolate but the original code was interpolating also for inac==2 so I kept it because I didn't want to mess up something.
Ok, thanks for the message and the comments. I will merge it (we need Baeck-An for a project here in Prague) and if you find any problems with implementation, we can fix it in other PR. |
|
Sounds good! Thanks for the final tweaks and documentation. |
Approximate Baeck-An couplings are calculated based on the history of PESs and allow to run FSSH without explicitly calculated couplings.