Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
29 changes: 29 additions & 0 deletions .github/workflows/fortitude.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
name: Lint with Fortitude

on:
workflow_dispatch:
pull_request:
branches: [main]
push:
branches: [main]

concurrency:
group: ${{ github.workflow }}-${{ github.event.pull_request.number }}
cancel-in-progress: true

jobs:
test:
name: Fortitude-Lint
runs-on: [ubuntu-latest]
strategy:
fail-fast: false
steps:
- uses: actions/checkout@v5
with:
fetch-depth: 0
- uses: actions/setup-python@v5
with:
python-version: '3.13'
- run: |
python -m pip install fortitude-lint
fortitude check --ignore C003 --output-format github
74 changes: 37 additions & 37 deletions Code/coulomb.f90
Original file line number Diff line number Diff line change
@@ -1,20 +1,20 @@
!------------------------------------------------------------------------------
! MODULE: Coulomb
!------------------------------------------------------------------------------
! DESCRIPTION:
! DESCRIPTION:
!> @brief
!!This module offers the subroutine \c poisson, which calculates the Coulomb
!!potential from the proton density. The two boundary conditions allowed are
!!distinguished by the value of the variable periodic. If it is true, the mesh
!!is assumed to be periodic in all three directions and the Poisson equation
!!is solved using the \f$\frac{1}{k^2}\f$ Green’s function in momentum space and
!!assuming a homogeneous negative background annulling the total charge (jellium
!!approximation). Otherwise the boundary condition is for an isolated charge
!!distribution. In this case the potential is calculated by the method of
!!doubling the grid in all directions with zero density, folding with the
!!This module offers the subroutine \c poisson, which calculates the Coulomb
!!potential from the proton density. The two boundary conditions allowed are
!!distinguished by the value of the variable periodic. If it is true, the mesh
!!is assumed to be periodic in all three directions and the Poisson equation
!!is solved using the \f$\frac{1}{k^2}\f$ Green’s function in momentum space and
!!assuming a homogeneous negative background annulling the total charge (jellium
!!approximation). Otherwise the boundary condition is for an isolated charge
!!distribution. In this case the potential is calculated by the method of
!!doubling the grid in all directions with zero density, folding with the
!!1/r-potential in momentum space and then restricting to the physical region.
!>
!>@details
!>@details
!!The Coulomb solver follows the ideas of \cite Eas79. We summarize
!!here the way it is realized in the code without detailed proofs.
!!
Expand All @@ -29,7 +29,7 @@
!!\\&
!! \mbox{coordinate spacing: } {\tt dx}, {\tt dy}, {\tt dz};
!!\\&
!! \mbox{momentum spacing: }
!! \mbox{momentum spacing: }
!! {\tt dk}_x
!! =
!! \frac{2\pi}{{\tt nx}\!\cdot\!{\tt dx}},
Expand All @@ -51,7 +51,7 @@
!! y_{\nu_y}=(\nu_y\!-\!{\tt ny}){\tt dy}\,,\,
!! z_{\nu_z}=(\nu_z\!-\!{\tt nz}){\tt dz}\,;
!!\\&
!! \mbox{momentum spacing: }
!! \mbox{momentum spacing: }
!! {\tt dk}_x
!! =
!! \frac{2\pi}{{2\tt nx}\!\cdot\!{\tt dx}},
Expand Down Expand Up @@ -80,19 +80,19 @@
!! directions this would lead to a value of \f$ 2.38/\Delta x \f$. Practical
!! experimentation, however, showed that the underlying assumption of a
!! point charge can be improved by some smeared-out density for nuclear
!! applications; a value of \f$ 2.84/\Delta x \f$ was found to be optimal.
!! We use the expression
!! applications; a value of \f$ 2.84/\Delta x \f$ was found to be optimal.
!! We use the expression
!! \f[ G(0) = \frac{2.84\sqrt{3}}{\sqrt{{\tt dx}^2+{\tt dy}^2+{\tt dz}^2}} \quad.\f]
!! -# Prepare the Greens function in \f$ k \f$-space by
!! 3D fast Fourier transformation (FFT) on the double grid \f$\mathcal{G}_2\f$.
!! \f$\tilde{G}(\mathbf{k}_{\bmu})={\tt FFT}\{G(\mathbf{r}_{\bnu})\}\f$.
!! The array \f$ \tilde{G}(\mathbf{k}_{\bmu}) \f$ is stored
!! The array \f$ \tilde{G}(\mathbf{k}_{\bmu}) \f$ is stored
!! for further continued use.
!!.
!!Once properly prepared, the practical steps for computing
!!the Coulomb field \f$ U_\mathrm{Coul}(\mathbf{r}_{\bnu}) \f$ for a
!!the Coulomb field \f$ U_\mathrm{Coul}(\mathbf{r}_{\bnu}) \f$ for a
!!density \f$\rho(\mathbf{r}_{\bnu})\f$ given on \f$\mathcal{G}_1\f$ are
!! -# Extend \f$ \rho_{\bnu} \f$ from \f$ \mathcal{G}_1 \f$ to \f$ \mathcal{G}_2 \f$
!! -# Extend \f$ \rho_{\bnu} \f$ from \f$ \mathcal{G}_1 \f$ to \f$ \mathcal{G}_2 \f$
!! by zero padding:
!! \f[ \rho_2(\mathbf{r}_{\bnu}) = \left\{\begin{array}{lcl} \rho(\mathbf{r}_{\bnu}) &\mbox{ if }& 1\leq \nu_i \leq {\tt n}i \mbox{ for } i\in\{x,y,z\}\\ 0 &\mbox{ else }\end{array}\right. \f]
!! -# Fourier transform the density in \f$ \mathcal{G}_2 \f$:
Expand All @@ -111,29 +111,29 @@ MODULE Coulomb
USE Params, ONLY: db,pi,e2
USE Grids, ONLY: nx,ny,nz,dx,dy,dz,wxyz,periodic
USE Densities, ONLY: rho
USE ISO_C_BINDING
USE ISO_C_BINDING, ONLY: c_long
IMPLICIT NONE
!>@name Dimensions of the grid on which the Fourier transform is calculated.
!>For the periodic case they are identical to the regular dimensions, for the
!>@name Dimensions of the grid on which the Fourier transform is calculated.
!>For the periodic case they are identical to the regular dimensions, for the
!>isolated case they are doubled.
!>@{
INTEGER,PRIVATE :: nx2,ny2,nz2
INTEGER,PRIVATE :: nx2,ny2,nz2
!>@}
!>@name Plans for FFTW complex forward and reverse transforms
!>@name Plans for FFTW complex forward and reverse transforms
!>with array dimensions depending on the boundary condition.
!>@{
INTEGER(C_LONG),PRIVATE,SAVE :: coulplan1,coulplan2
!>@}
!>@}
REAL(db),ALLOCATABLE,SAVE :: wcoul(:,:,:) !< the Coulomb potential as a three-dimensional array. Units: MeV.
COMPLEX(db),PRIVATE,ALLOCATABLE,SAVE :: q(:,:,:) !<array for the complex Green’s function (isolated) or array
!!of 1/r values. Its dimension also depends on the boundary condition.
PUBLIC :: poisson,coulinit,wcoul
PRIVATE :: initiq
PRIVATE
CONTAINS
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
! DESCRIPTION: poisson
!> @brief
!!This subroutine solves the Poisson equation by the Fourier method.
!!This subroutine solves the Poisson equation by the Fourier method.
!>
!> @details
!!For the periodic case, this means to just Fourier transform, multiply by
Expand All @@ -147,7 +147,7 @@ MODULE Coulomb
!!The key to success is not to use simply \f$1/k^2\f$ for the Fourier
!!transform of \f$1/r\f$, but to compute it on the actual grid the same way
!!as the densities and potentials are transformed
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
SUBROUTINE poisson(rhocharge)
REAL(db),INTENT(IN) :: rhocharge(nx,ny,nz)
COMPLEX(db),ALLOCATABLE :: rho2(:,:,:)
Expand All @@ -169,25 +169,25 @@ SUBROUTINE poisson(rhocharge)
wcoul=REAL(rho2(1:nx,1:ny,1:nz))/(nx2*ny2*nz2)
DEALLOCATE(rho2)
END SUBROUTINE poisson
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
! DESCRIPTION: coulinint
!> @brief
!!This subroutine does the necessary initialization.
!!This subroutine does the necessary initialization.
!!
!>
!> @details
!!It calculates the dimension for the Fourier transform depending on
!!the boundary condition, allocates the necessary arrays and sets up
!!It calculates the dimension for the Fourier transform depending on
!!the boundary condition, allocates the necessary arrays and sets up
!!the FFTW plans.
!!Then it composes the array \c q. For the periodic case this
!!contains the values of the Green's function \f$ 1/k^2 \f$, with zero at the
!!origin, while for the isolated case it first calculates the inverse
!!shortest distance $1/r$ from the origin with indices (1,1,1), replaces
!!the value at the origin and then
!!Fourier-transforms this to use it for folding the density with the
!!coordinate-space Green's function.
!!coordinate-space Green's function.
!>
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
SUBROUTINE coulinit
INCLUDE 'fftw3.f'
REAL(db),ALLOCATABLE :: iqx(:),iqy(:),iqz(:)
Expand Down Expand Up @@ -223,11 +223,11 @@ SUBROUTINE coulinit
END IF
DEALLOCATE(iqx,iqy,iqz)
END SUBROUTINE coulinit
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
! DESCRIPTION: initiq
!> @brief
!!This subroutine calculates the contributions of each Cartesian index
!!to \f$r^2\f$ and \f$k^{-2}\f$.
!!to \f$r^2\f$ and \f$k^{-2}\f$.
!>
!> @details
!!This depends on the boundary condition. For
Expand All @@ -248,7 +248,7 @@ END SUBROUTINE coulinit
!> REAL(db), takes the grid spacing.
!> @param[out] iq
!> REAL(db), array, returns the values for \f$r^2\f$ and \f$k^{-2}\f$.
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
SUBROUTINE initiq(n,d,iq)
INTEGER,INTENT(IN) :: n
REAL(db),INTENT(IN) :: d
Expand Down
38 changes: 20 additions & 18 deletions Code/densities.f90
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
!------------------------------------------------------------------------------
! MODULE: Densities
!------------------------------------------------------------------------------
! DESCRIPTION:
! DESCRIPTION:
!> @brief
!!This module has two purposes: it defines and allocates the densities
!!and currents making up the mean field,
Expand All @@ -16,8 +16,10 @@ MODULE Densities
USE Levels, ONLY: cdervx,cdervy,cdervz
USE Trivial, ONLY: cmulx,cmuly,cmulz
IMPLICIT NONE
PRIVATE
PUBLIC :: rho,tau,sdens,sodens,current,alloc_densities
SAVE
!>@name Scalar densities:
!>@name Scalar densities:
!>These are dimensioned <tt>(nx,ny,nz,2)</tt>,
!>where the last index is 1 for neutrons and 2 for protons.
!>@{
Expand All @@ -34,7 +36,7 @@ MODULE Densities
!!\f$ \hbar^2/2m \f$. Units: \f$ {\rm fm}^{-5}$.
!>@}
!
!>@name Vector densities:
!>@name Vector densities:
!>These are dimensioned <tt>(nx,ny,nz,3,2)</tt>, where the last index is 1 for neutrons and 2 for
!>protons, and the next-to-last stands for the Cartesian direction.
!>@{
Expand All @@ -50,7 +52,7 @@ MODULE Densities
REAL(db),ALLOCATABLE,DIMENSION(:,:,:,:,:) :: sdens
!<the spin density. It is defined as
!!\f[ \vec\sigma_q(\vec r)=\sum_{\alpha\in q} w_\alpha^2 \sum_{ss'}\psi_\alpha^*(\vec
!!r,s)\,\sigma_{ss'} \,\psi_\alpha(\vec r,s').\f]
!!r,s)\,\sigma_{ss'} \,\psi_\alpha(\vec r,s').\f]
!!Note that it does not include the factor \f$ \hbar/2 \f$. Units: \f$ {\rm fm}^{-3} \f$.
REAL(db),ALLOCATABLE,DIMENSION(:,:,:,:,:) :: sodens
!<the spin-orbit density, defined as
Expand All @@ -60,16 +62,16 @@ MODULE Densities
!!Its units are also \f$ {\rm fm}^{-4} \f$.
!>@}
CONTAINS
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
! DESCRIPTION: alloc_densities
!> @brief
!!This is simply a short routine to allocate all the arrays defined in this module.
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
SUBROUTINE alloc_densities
ALLOCATE(rho(nx,ny,nz,2),tau(nx,ny,nz,2),current(nx,ny,nz,3,2), &
sdens(nx,ny,nz,3,2),sodens(nx,ny,nz,3,2))
END SUBROUTINE alloc_densities
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
! DESCRIPTION: add_density
!> @brief
!!This subroutine is given a single-particle wave function \c psin
Expand All @@ -93,7 +95,7 @@ END SUBROUTINE alloc_densities
!!The separate copies are then combined using the \c OPENMP <tt>REDUCE(+)</tt> directive.
!!
!!The local copies of the densities passed as arrays are denoted with
!!the prefixed letter "l" for \a local; they are \c lrho, \c ltau,
!!the prefixed letter "l" for \a local; they are \c lrho, \c ltau,
!!\c lcurrent, \c lsdens, and \c lsodens.
!!
!!If the weight is zero, there is nothing to do and the subroutine
Expand All @@ -110,7 +112,7 @@ END SUBROUTINE alloc_densities
!!
!!The complex products always in the end evaluate to something real and
!!the expressions are simplified to take this into account. For example,
!!the following transformation is done:
!!the following transformation is done:
!!\f{eqnarray*}{
!!\frac{1}{2\I}(\psi^*\nabla\psi-\psi\nabla\psi^*)&=&
!!\frac{1}{2\I}\left(\psi^*\nabla\psi-(\psi^*\nabla\psi)^*\right)\\
Expand Down Expand Up @@ -140,14 +142,14 @@ END SUBROUTINE alloc_densities
!> REAL(db), array, takes and adds the spin density.
!> @param[in, out] lsodens
!> REAL(db), array, takes and adds the spin-orbit density.
!---------------------------------------------------------------------------
SUBROUTINE add_density(iq,weight,psin,lrho,ltau,lcurrent,lsdens,lsodens)
!---------------------------------------------------------------------------
SUBROUTINE add_density(iq,weight,psin,lrho,ltau,lcurrent,lsdens,lsodens)
COMPLEX(db),INTENT(INOUT) :: psin(nx,ny,nz,2)
REAL(db),DIMENSION(:,:,:,:),INTENT(INOUT) :: lrho,ltau
REAL(db),DIMENSION(:,:,:,:,:),INTENT(INOUT) :: lcurrent,lsdens,lsodens
INTEGER,INTENT(IN) :: iq
REAL(db),INTENT(IN) :: weight
COMPLEX(db),ALLOCATABLE :: ps1(:,:,:,:)
COMPLEX(db),ALLOCATABLE :: ps1(:,:,:,:)
INTEGER :: ix,iy,iz
ALLOCATE(ps1(nx,ny,nz,2))
IF(weight<=0.D0) RETURN
Expand All @@ -170,9 +172,9 @@ SUBROUTINE add_density(iq,weight,psin,lrho,ltau,lcurrent,lsdens,lsodens)
! x-derivatives
!***********************************************************************
IF(TFFT) THEN
CALL cdervx(psin,ps1)
CALL cdervx(psin,ps1)
ELSE
CALL cmulx(der1x,psin,ps1,0)
CALL cmulx(der1x,psin,ps1,0)
ENDIF
ltau(:,:,:,iq)=ltau(:,:,:,iq)+weight* &
(ps1(:,:,:,1)*CONJG(ps1(:,:,:,1))+ps1(:,:,:,2)*CONJG(ps1(:,:,:,2)))
Expand All @@ -191,9 +193,9 @@ SUBROUTINE add_density(iq,weight,psin,lrho,ltau,lcurrent,lsdens,lsodens)
! y-derivatives
!***********************************************************************
IF(TFFT) THEN
CALL cdervy(psin,ps1)
CALL cdervy(psin,ps1)
ELSE
CALL cmuly(der1y,psin,ps1,0)
CALL cmuly(der1y,psin,ps1,0)
ENDIF
ltau(:,:,:,iq)=ltau(:,:,:,iq)+weight* &
(ps1(:,:,:,1)*CONJG(ps1(:,:,:,1))+ps1(:,:,:,2)*CONJG(ps1(:,:,:,2)))
Expand All @@ -212,9 +214,9 @@ SUBROUTINE add_density(iq,weight,psin,lrho,ltau,lcurrent,lsdens,lsodens)
! z-derivatives
!***********************************************************************
IF(TFFT) THEN
CALL cdervz(psin,ps1)
CALL cdervz(psin,ps1)
ELSE
CALL cmulz(der1z,psin,ps1,0)
CALL cmulz(der1z,psin,ps1,0)
ENDIF
ltau(:,:,:,iq)=ltau(:,:,:,iq)+weight* &
(ps1(:,:,:,1)*CONJG(ps1(:,:,:,1))+ps1(:,:,:,2)*CONJG(ps1(:,:,:,2)))
Expand Down
Loading