Good First Issues in Fortran

Brief summary of bug There are two FATES boundary conditions relating to liquid water that that are correctly set in wrap_btran() but set to total soil water in dynamics_driv(). The incorrect dynamics_driv() calls overwrite the correct wrap_btran() calls, causing incorrect boundary conditions, but only in the first timestep of the day—see where those subroutines are called relative to each other in clm_drv(). (Full disclosure: This bug was identified by Claude during my moss work, but I've traced it myself to check. I've also typed this issue out myself.) General bug information CTSM version you are using: ctsm5.4.054 Does this bug cause significantly incorrect results in the model's science? Probably not meaningfully so, but answers are affected. Configurations affected: Science affected in FATES runs with soil ice present: Runs with any PFT that is "stress (hydro) deciduous" (4 in default FATES parameter file) or "stress (hydro) semi-deciduous" (none). Maybe others, but that first

Should we maybe call this something like time_sampling_threshold or something like that? I wonder if there could be some confusion with horiz_remap_fill_threshold. I'm not opposed to the current name, just checking if there are concerns about confusion. Originally posted by @bartgol in #8481 (comment)

Dear cp2k developers, I have noticed some reported works but have fallen into the inappropriate descriptions about the use of constrained AIMD (cAIMD) simulation for the achievement of PMF. As the latest revised description for the cAIMD, the blue moon ensemble average connects constrained and unconstrained molecular dynamics beyond the Lagrange multiplier associated with the parameter used in the SHAKE algorithm rather than the pure Lagrange multiplier. Here, these works directly interpret the Lagrange multiplier (λ) from SHAKE algorithm as the ensemble-averaged constrained force, while absolutely neglects the correction from Z and G. In these works, this enhanced sampling is used to obtain the activation energy, but I think sometimes this is very dangerous to crudely neglect the correction from Z and G beyond λ. The CP2K's official example is OK (https://www.cp2k.org/exercises:2015_ethz_mmm:nacl_free_energy), because for the example of NaCl dissociation, the CV is only the simple d

Dear developers, It seems that some questions for the calculation of hirshfeld in cp2k are available: https://groups.google.com/g/cp2k/c/bf0RsWg5Ty8/m/J6y6-EDpAwAJ https://groups.google.com/g/cp2k/c/Vo6kuq6UCnc/m/rEVkjqZcCwAJ The question lies in that the method seems not the follow of standard scheme. The SHAPE_FUNCTION GAUSSIAN frequently output the charges that does not follow the normal chemical meaning. And for the SHAPE_FUNCTION DENSITY with SELF_CONSISTENT T cannot be strictly attributed into HIRSHFELD-I strategy. If any developers are resolving this confusion? Maybe we can try to retain the old definition in the code and then add some tags for newly definitions.

CP2K currently has several Libint-related layers with overlapping and implicit responsibilities. This may make the existing code difficult to extend when adding new electronic-integral APIs or alternative libraries. At present, CP2K's Libint integration is distributed across several files, including src/libint_wrapper.F, src/libint_2c_3c.F and src/hfx_libint_interface.F. I am wondering whether a more unified Libint interface would be useful and, if so, how its responsibilities might be organized. Related Issues: #5741 #1232

CP2K currently contains several recurrence implementations for Gaussian integrals, including: coulomb2 in src/aobasis/ai_coulomb.F; coulomb2_new in src/aobasis/ai_coulomb.F; coulomb3 in src/aobasis/ai_coulomb.F; and many other APIs. At the same time, CP2K already provides libint wrappers, which raised a question: should some of the existing APIs migrated to an external integral backend, such as libint or libcint mentioned in #1232?

The scalapack.pc installed by tools/toolchain (ScaLAPACK 2.2.3) has libdir=/lib instead of ${prefix}/lib, and Requires: mpi lapack blas. The latter makes pkg-config resolve BLAS/LAPACK transitively against the system .pc files, ignoring the OpenBLAS the toolchain just built and whatever -DCP2K_BLAS_VENDOR was set to. ssmp builds fine; only the MPI versions pull in ScaLAPACK. Environment CP2K version/commit: d734930 Debian 13 (trixie), gcc 14.2.0, system OpenMPI, toolchain-installed openblas-0.3.34 + scalapack-2.2.3 ./install_cp2k_toolchain.sh -j 16 --with-plumed=install --with-deepmd=install --with-libtorch=install Installed file tools/toolchain/install/scalapack-2.2.3/lib/pkgconfig/scalapack.pc: prefix=/path/to/cp2k/tools/toolchain/install/scalapack-2.2.3 libdir=/lib Name: scalapack Version: 2.2.3 Libs: -L -lscalapack Requires: mpi lapack blas Note libdir=/lib and the bare -L. Result $ pkg-config --libs lapack blas -L/usr/lib/x86_64-linux-gnu/openblas-openmp -llapack -lblas FindSCA

Euler–Lagrange two-way coupling gives rank-count-dependent results — investigated: round-off seed amplified by chaos, not a decomposition bug Summary Moving Lagrangian bubbles with two-way coupling (lag_params%vel_model = 2, lag_params%solver_approach = 2) produce visibly different results when the same case is run with different MPI rank counts. A systematic bisection campaign (attached tarball) traced this to a one-ulp (2.2e-16) summation-order difference in the two-way smearing deposition wherever a bubble's kernel straddles a subdomain boundary — the halo-accumulate reduce in s_mpi_reduce_beta_variables_buffers sums the same contributions in a decomposition-dependent order. The strongly chaotic coupled-bubble dynamics then amplify that seed exponentially (measured e-folding time ≈ 0.38 nondim time units) to O(1) trajectory differences within a few hundred steps. Verdict: there is no decomposition bug in the moving-bubble / two-way-coupling path. Every mechanism tested in is

Summary s_compute_viscous_stress_tensor (added in #1614 for higher-order central-difference IB drag) performs an out-of-bounds heap read whenever an immersed-boundary body's finite-difference stencil is centered on a cell adjacent to a domain boundary. The bug has been latent in master since #1614 (2026-06-24); it is silent in optimized builds and only manifests as a hard crash under bounds-checking (./mfc.sh ... --debug/reldebug) builds. cc @danieljvickers (author of #1614 and #1618) Root cause In src/simulation/m_viscous.fpp, the velocity-gradient stencil indexes the FD-coefficient arrays by cell: velocity_gradient_tensor(l, 2) = velocity_gradient_tensor(l, 2) & + fd_coeff_y(r, j)*q_prim_vf(...)%sf(i, j + r, k) But fd_coeff_x/y/z are allocated over the interior only — e.g. fd_coeff_y(-fd_number:fd_number, 0:n) (in m_derived_variables.fpp). The IB drag path (s_compute_ib_forces in m_ibm.fpp) evaluates this stencil at body ghost-point cells, which for a body sitting against a bou

Summary Add a helper (toolchain utility and/or smarter case-file parameterization) to make hyperbolic-tangent grid-stretching parameters easier to choose. Today users have to hand-tune the raw stretch parameters (a_x, x_a, x_b, loops_x, etc.) by trial and error, which is unintuitive and error-prone. Motivation This came out of the discussion on #1625, which fixed MPI-rank-dependent divergent grid stretching. Now that stretching is rank-invariant, it's a good time to make setup friendlier. As @danieljvickers suggested: I think it may also be neat to add a grid-stretching helper to the toolchain or a smarter parameter grid stretching in the case files to make getting your parameters more obvious. Something to help people with parameter selection based on their desired domain size, uniform region, uniform spacing, etc. Proposed idea Let users specify what they actually care about, and derive the raw stretch parameters from it: desired total domain size extent of the uniform (unstretche

In the recent addition of higher-order finite differencing for immersed boundaries from #1614, the fd_order case file parameter has begun to be overloaded in terms of its use. 4 different case file settings all use the same finite difference order. This means that we are somewhat less flexible and the term appears more ambiguous to users. It would be helpful to consider a refactor of the finite difference scheme for paramters like probes, hypo-elasticity, immersed boundaries, MHD, post-processing, etc. In particular, the finite difference order used in all of these terms should be disambiguated, The goal seems to be separate the use cases of this in pre_process, simulation, post_process to be distinct

Brief summary of bug $ grep -R -n qflx_evap_can_col ./* ./biogeophys/WaterFluxType.F90:48: real(r8), pointer :: qflx_evap_can_col (:) ! col evaporation from leaves and stems (mm H2O/s) (+ = to atm) ./biogeophys/WaterFluxType.F90:239: call AllocateVar1d(var = this%qflx_evap_can_col, name = 'qflx_evap_can_col', & General bug information CTSM version you are using: master/latest Does this bug cause significantly incorrect results in the model's science? No Configurations affected: All Details of bug This should not even be allocated ever. There is qflx_evap_veg_col Important details of your setup / configuration so we can reproduce the bug Important output or errors that show the problem

The IB patches have just been refactored to significantly reduce the number of lines of code required in the module. This was done by making use of shared subroutines and functions that are common to all IB patch set up, with mild if/then trees for various patches. There is arguably more overlapping/redundant code in the ICPP patches than there was in the IB patch code, and similar improvements can be made. The product of this should be unifying as many subroutines as possible to use common initial set up and geometry checking. The result should be a significant reduction in the number of lines of code in the module.

The IB patches have a rotation system that has been fully integrated across all patches, and the ICPP patch code now makes use of underlying functions compatible with rotation. MFC would be enhanced by integrating rotations into the ICPP patch geometries greatly expanding the flexibility of case generation without requiring hard-coded patches. This feature should Add a matching rotation framework to icpp patches to be the same as the patch_ib rotation code Move the rotation framework to the m_patch_geometries module add optimizations to allow faster projection of icpp patches onto the grid, such as the bounding indices used in IB patches. At least one example of a rotated patch

We use various integers to differentite between geometries in 2D and 3D. For example, 2D STL models and 3D STL models are differentiated. Airfoils, and more are also included. Use of the code may be simpler if these patch geometry definitions were unified into a single value. I propose the following use the same patch geometry values: 2D and 3D STL models 2D and 3D airfoils circles and cylinders (or spheres depending on discretion) rectangles and cuboids To maintain legacy functionality, we can have the 3D cases support either value.

When wee run test with the --only flag, it is case sensitive. This is slightly frustrating because flags include "STL", "IBM", "Sphere", "Cuboid", "slip", etc. Since there is no convention in the way test labels are captialized, it requires the used to remember the specific capitalization for the test that they want to run. It would make much more sense to not be case-sensitive. This issue can be solved by modifications to the toolchain for the ./mfc.sh test --only <label> feature. When complete, the test suite should run the same with various spellings: STL, Stl, stl, sTl`, etc.

Problem it solves When one runs the precheck command in the toolchain, one of the steps runs validate on all examples and reports the number of examples that fail, but not which examples. This requires the used to check cases one-at-a-time manually. One can accomplish this with BASH commands, but it would be much nicer if the validator had something like an --all flag and just tries to validate all examples. Proposed solution Add an --all flag to validator that checks all examples. The flag should accurately report which cases fail validation the precheck for commits should swap from what it uses now and instead use this command (Optional) Validation is also very slow, which leads me to always commit with --no-verify. It would be nice if it were much faster. Alternatives You could also simply let the validate command accept wildcard values, like *, or directories and the user could simply pass in something like ./mfc.sh validate examples or ./mfc.sh validate examples/* to accomplish

I have done some test runs in single precision for the ABL case (using examples/shear_convection_abl and other similar setups). It seems to be broken. This issue is to track this problem and hopefully solve it. Fluid stats: Temperature stats: I don't know exactly how to tackle this, so here is some random observations: the variable that really seems to be broken is temperature. I am thinking that everything "starts" there and the propagates to the velocity fields via the Boussinesq term in the momentum equations (in ABL cases, that's how temperature is coupled to velocity) however, there are some weird stuff in the fluid statistics as well: see for instance the little wiggles above 750m in it doesn't seem to be dependent on resolution. I have tested both a coarser and finer resolution and the problem is the same the equivalent Double Precision case works perfectly the tests have been carried out on LUMI only I am thinking to start debugging this by running without temperature firs

Several examples (and tests) make use of SimpleNet, and we also have MultiIONet and BatchingNet. To reduce repetition it might be useful to move these into a central part of ftorch-utils as single points of truth. Original discussion in #555 (comment)

From my perspective, the current CP2K documentation could be further improved in terms of usability, especially for new users. A key issue is that many keywords in the Input Reference are described only very briefly, without sufficient explanation of their practical meaning and usage in real calculations. Important aspects such as typical use cases, parameter sensitivity, and common pitfalls are often not covered. As a result, users frequently need to consult the original literature to understand the underlying methods and implementation details, and rely on repeated trial and error to use certain features effectively. This leads to a fairly steep learning curve for CP2K. I'm aware that the current Input Reference is directly generated from source code, and one effective approach would be to extend these descriptions in the Fortran source. However, this may not always be the most suitable solution, especially when more detailed, tutorial-style explanations are needed. As an alternative

... like min, max, avg, std, var, etc. also mahf708:hbc2 has some expanded capabilities in vert reductions diag notes: this can be done as a rigid "specify, then spit out" like T_mid_temp_avg_5_[TIME_FREQ] and then delegate to IO to output as usual or better, it can be done as a diag for TIME-RUNNING (or rolling) statistics, which is a good anticipation for the streaming future, like T_mid_temp_avg always saves the average and IO can request it whenever it wants or even better, it can be done as a mix of both by inserting knowledge of the last IO action (from the output manager/controller) into the diag itself (since the io manager spawns out a specific diag instance)

This is a collection of ideas and an extension of an older threads #225 , #4811 OT SCCS DFT+U (Mulliken, Mulliken charges, Lowdin) DFT+U (Tensorial) PDOS k-point symmetry reduction Wannier functions, interface to Wannier90 TDDFT alternative diagonalization schemes for k-points: LANCZOS, FILTER_MATRIX #496 restart k-point calculation from Gamma-point calculation Alternative HFX (for instance two-center RI as used in FHI-aims or (partially) reciprocal-space scheme as used in PySCF, Head-Gordon group) post-SCF-methods (MP2, RPA, GW) Extrapolation schemes of wavefunction in MD or opt jobs, like PS and ASPC as specified by &FORCE_EVAL/&DFT/&QS/EXTRAPOLATION Periodic electric field as requested by &FORCE_EVAL/&DFT/&PERIODIC_EFIELD (which may or may not be a separate task from OT) Linear response properties, like NMR, as specified by &FORCE_EVAL/&PROPERTIES/&LINRES Löwdin/Loewdin population analysis as requested by &FORCE_EVAL/&DFT/&PRINT/&LOWDIN cp2k/src/pop

After PR #4799, CP2K can now be compiled successfully with Intel oneAPI 2025.2.1 (base Toolkit with MKL, and HPC Toolkit with icx, icpx, ifx and MPI library). However, during a test of cell optimization task, I found the Intel-compiled CP2K runs not as fast as what's compiled using GCC with OpenMPI or MPICH. Even more strangely, there occurs a memory issue as shown below. (increasing slowly and abnormally) When I try falling Intel oneAPI version back to 2025.2.1, the memory leak still occurs. While here's the condition with GCC/OpenMPI, which is normal and uses only half of the memory compared to the Intel-compiled one: It's also worth noticing that regtests with ifx are currently disabled with jobs requiring libint for runtime failure and wrong results, and optimize level of several files has to be reduced. Therefore, in addition to compilation-related issues, I would like to further discuss whether there might be some compatibility issues between the current CP2K code and oneAPI.

In #489 we moved to rename the test utils assert_isclose and assert_allclose as isclose and allclose. However, we maintained aliases for the original versions for backwards compatibility, with a deprecation warning.

Because this is confusing: CTSM/src/biogeochem/CNPhenologyMod.F90 Line 2095 in 19d9b8a leafout => crop_inst%gddtsoi_patch , & ! Input: [real(r8) (:) ] gdd from top soil layer temperature Seemingly this is the only place it's used.

We commonly need to check that soil layers are oversaturated and, if so, move the excess soil water up by one soil layer (or runoff if the top layer is oversaturated). This happens in least three places. This should be modularized to avoid repeated code.

There is a variety of features without any coverage by the regression tests (either files or functionalities from the input file). Have they ever been practically used in the past or do you know whether they are still functionable? From the older features, it is (incomplete list) RIGPW ET_COUPLING%TYPE_OF_CONSTRAINT DDAPC src/ct_methods.F (Cayley transforms) src/kpsym.F (Tools for k-point symmetry) src/minbas_methods.F NEGF%INTEGRATION_METHOD CLENSHAW-CURTIS, NEGF EMBED%EMBED_METHOD DMFET PEXSI DFT%SCF%OUTER_SCF%TYPE [S2_CONSTRAINT|BASIS_CENTER_OPT] DF%PRINT%AO_MATRIX%COMMUTATOR_HR VCD_env_type%do_mpf (not set from input anywhere) DFT%PRINT%WANNIER90 se_fock_matrix_dbg.F (never called by caller) QS%SE%PERIODIC EWALD_R3 (there is only an untested regtest) src/motion/free_energy_methods.F: ui_check_convergence (test to call it never passes) QS/regtest-tipscan (was not registered) DFT%TRANSPORT DFT%PRINT%WANNIER90 INTERPOLATOR%KIND [SPLINE3|SPLINE3_NOPBC] (does definit