From 570b3b804b58557b477ecdcfc427035d095a09b3 Mon Sep 17 00:00:00 2001 From: mcbarton Date: Wed, 17 Sep 2025 15:16:38 +0100 Subject: [PATCH 01/17] Add Fortitude linter workflow (#2) --- .github/workflows/fortitude.yml | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 .github/workflows/fortitude.yml diff --git a/.github/workflows/fortitude.yml b/.github/workflows/fortitude.yml new file mode 100644 index 0000000..15f0c3d --- /dev/null +++ b/.github/workflows/fortitude.yml @@ -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 --output-format github From 15f4a432928add4f157b0e7757ed23ebe82f42a9 Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Wed, 17 Sep 2025 15:19:25 +0100 Subject: [PATCH 02/17] Fix Fortitude rule S101 using automatic fix flag --- Code/coulomb.f90 | 70 ++++---- Code/densities.f90 | 36 ++-- Code/dynamic.f90 | 132 +++++++------- Code/energies.f90 | 32 ++-- Code/external.f90 | 116 ++++++------ Code/forces.f90 | 40 ++--- Code/fourier.f90 | 8 +- Code/fragments.f90 | 112 ++++++------ Code/grids.f90 | 90 +++++----- Code/inout.f90 | 106 +++++------ Code/levels.f90 | 74 ++++---- Code/main3d.f90 | 8 +- Code/meanfield.f90 | 138 +++++++------- Code/moment.f90 | 140 +++++++-------- Code/pairs.f90 | 244 ++++++++++++------------- Code/parallel.f90 | 46 ++--- Code/params.f90 | 20 +-- Code/sequential.f90 | 10 +- Code/static.f90 | 280 ++++++++++++++--------------- Code/trivial.f90 | 132 +++++++------- Code/twobody.f90 | 330 +++++++++++++++++----------------- Code/user.f90 | 18 +- Code/ylm.f90 | 22 +-- Test/GR/spectral_analysis.f90 | 12 +- Utils/Cuts/Cuts.f90 | 2 +- Utils/Fileinfo/fileinfo.f90 | 2 +- Utils/Overlap/overlap.f90 | 2 +- 27 files changed, 1111 insertions(+), 1111 deletions(-) diff --git a/Code/coulomb.f90 b/Code/coulomb.f90 index 5c111cc..a06bbc7 100644 --- a/Code/coulomb.f90 +++ b/Code/coulomb.f90 @@ -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. !! @@ -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}}, @@ -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}}, @@ -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$: @@ -113,27 +113,27 @@ MODULE Coulomb USE Densities, ONLY: rho USE ISO_C_BINDING 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(:,:,:) ! @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 @@ -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(:,:,:) @@ -169,15 +169,15 @@ 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 @@ -185,9 +185,9 @@ END SUBROUTINE poisson !!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(:) @@ -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 @@ -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 diff --git a/Code/densities.f90 b/Code/densities.f90 index 1b9482d..7de24db 100644 --- a/Code/densities.f90 +++ b/Code/densities.f90 @@ -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, @@ -17,7 +17,7 @@ MODULE Densities USE Trivial, ONLY: cmulx,cmuly,cmulz IMPLICIT NONE SAVE - !>@name Scalar densities: + !>@name Scalar densities: !>These are dimensioned (nx,ny,nz,2), !>where the last index is 1 for neutrons and 2 for protons. !>@{ @@ -34,7 +34,7 @@ MODULE Densities !!\f$ \hbar^2/2m \f$. Units: \f$ {\rm fm}^{-5}$. !>@} ! - !>@name Vector densities: + !>@name Vector densities: !>These are dimensioned (nx,ny,nz,3,2), where the last index is 1 for neutrons and 2 for !>protons, and the next-to-last stands for the Cartesian direction. !>@{ @@ -50,7 +50,7 @@ MODULE Densities REAL(db),ALLOCATABLE,DIMENSION(:,:,:,:,:) :: sdens !@} 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 @@ -93,7 +93,7 @@ END SUBROUTINE alloc_densities !!The separate copies are then combined using the \c OPENMP REDUCE(+) 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 @@ -110,7 +110,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)\\ @@ -140,14 +140,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 @@ -170,9 +170,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))) @@ -191,9 +191,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))) @@ -212,9 +212,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))) diff --git a/Code/dynamic.f90 b/Code/dynamic.f90 index 80c3d80..73c6e71 100644 --- a/Code/dynamic.f90 +++ b/Code/dynamic.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! MODULE: Dynamic !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module contains the routines needed for time propagation of the !!system. All the logic needed for this @@ -31,26 +31,26 @@ MODULE DYNAMIC !! expansion of the potential INTEGER :: mrescm=0 !< frequency of c.m. motion correction REAL(db) :: rsep !< the final separation distance. The calculation - !! is stopped if there has been a reseparation into two fragments and their distance exceeds \c rsep. + !! is stopped if there has been a reseparation into two fragments and their distance exceeds \c rsep. LOGICAL :: texternal=.FALSE. !< this logical variable indicates that an !! external field is present. See module \c External. LOGICAL :: text_timedep !< this logical variable indicates that the !! external field is time-dependent and does not describe an instantaneous boost. - REAL(db),PARAMETER :: esf=0.0D0 !< this is the energy shift for the call to \c hpsi. + REAL(db),PARAMETER :: esf=0.0D0 !< this is the energy shift for the call to \c hpsi. !! Since it is not used in the dynamics part of the code, it !! is here set to the constant value of zero. CONTAINS -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: getin_dynamic !> @brief !!This is a relatively simple routine that reads the input for namelist !!\c dynamic and prints it on standard output. If \c texternal is !!true, it also calls \c getin_external to read the external field !!parameters. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE getin_dynamic NAMELIST /dynamic/ nt,dt,mxpact,mrescm,rsep,texternal - READ(5,dynamic) + READ(5,dynamic) IF(wflag) THEN WRITE(*,*) '***** Parameters for the dynamic calculation *****' WRITE(*,"(A,I8,A,F10.5,A)") " Number of time steps:",nt, & @@ -61,7 +61,7 @@ SUBROUTINE getin_dynamic ENDIF IF(texternal) CALL getin_external END SUBROUTINE getin_dynamic -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: dynamichf !> @brief !!This subroutine performs the main time-integration algorithm starting @@ -69,25 +69,25 @@ END SUBROUTINE getin_dynamic !> !> @details !!Its building blocks are: -!! - Step 1: preparation phase: this phase consists of several substeps. +!! - Step 1: preparation phase: this phase consists of several substeps. !! -# If this is not a restart, the time and iteration number are !! zeroed (for a restart only the physical time is taken from the !! \c wffile). The wave functions are saved in the !! \c wffile, to save setup time in case the calculation has to !! be restarted from this initial point. -!! -# The instantaneous external boost is applied using -!! \c extboost. If this subroutine has applied a boost, it sets +!! -# The instantaneous external boost is applied using +!! \c extboost. If this subroutine has applied a boost, it sets !! \c text_timedep to \c .FALSE. so no further calls to !! external-field routines are made (except for \c print_extfield). !! -# The protocol files *.res are initialized with -!! their header lines. +!! their header lines. !! -# The densities and current are calculated in a loop over wave -!! function by calling \c add_densities. +!! function by calling \c add_densities. !! They are first set to zero and then accumulated in a !! loop over the set on the local node, followed by collecting them !! over all nodes. !! -# The mean field and (if \c texternal is true) the -!! external field are calculated for time zero using routines +!! external field are calculated for time zero using routines !! \c skyrme and \c extfld. !! -# Then \c tinfo is called to calculate and print the !! single-particle quantities and the total energies at time 0 or @@ -100,7 +100,7 @@ END SUBROUTINE getin_dynamic !! !! - Step 2: predictor time step: the loop over the iteration !! index \c iter is started. Then the densities and mean-field -!! components are estimated, i. e., in effect the Hamiltonian +!! components are estimated, i. e., in effect the Hamiltonian !! \f$ \hat h(t+\tfrac1{2}\Delta t) \f$ required by the numerical method. This is done !! by evolving the wave functions for a full \c dt using the old !! Hamiltonian and averaging the densities between old and new ones to @@ -135,16 +135,16 @@ END SUBROUTINE getin_dynamic !! -# The densities are reset to zero before the wave function loop, !! so the densities summed up are the purely the densities at the end !! of the time step, -!! -# the series expansion in \c tstep now uses the full +!! -# the series expansion in \c tstep now uses the full !! \c mxpact terms, and !! -# the new wave functions are copied back into \c psi to be !! available for the next time step. !! !! - Step 4: Center-of-mass correction: !! If a center-of-mass correction is desired by the user by setting -!! mrescm /=0, +!! mrescm /=0, !! subroutine \c resetcm is called every \c mrescm'th time step to -!! reset the center-of-mass velocity to zero. +!! reset the center-of-mass velocity to zero. !! !! - Step 5: generating some output At this point the time is !! advanced by \c dt because the physical time is now the end of @@ -159,7 +159,7 @@ END SUBROUTINE getin_dynamic !! onto \c wffile depending on \c mrest. !! !!This ends the time loop and subroutine \c dynamichf itself. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE dynamichf INTEGER :: nst,istart COMPLEX(db),ALLOCATABLE :: ps4(:,:,:,:) @@ -167,7 +167,7 @@ SUBROUTINE dynamichf ! Step 1: Preparation phase IF(.NOT.trestart) THEN iter=0 - time=0.0D0 + time=0.0D0 ! save wave functions CALL write_wavefunctions IF(wflag) WRITE(*,*) 'Wrote wave function file after initialization' @@ -216,7 +216,7 @@ SUBROUTINE dynamichf !$OMP REDUCTION(+:rho,tau,current,sdens,sodens) DO nst=1,nstloc CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), & - psi(:,:,:,:,nst),rho,tau,current,sdens,sodens) + psi(:,:,:,:,nst),rho,tau,current,sdens,sodens) ENDDO !$OMP END PARALLEL DO IF(tmpi) CALL collect_densities @@ -227,7 +227,7 @@ SUBROUTINE dynamichf !*********************************************************************** istart=iter+1 ! Step 2: start loop and do half-time step - Timestepping: DO iter=istart,nt + Timestepping: DO iter=istart,nt IF(wflag) WRITE(*,'(/A,I6,A,F8.2,A)') ' Starting time step #',iter, & ' at time=',time,' fm/c' ! correction for parallel version @@ -242,10 +242,10 @@ SUBROUTINE dynamichf !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,ps4) SCHEDULE(STATIC) & !$OMP REDUCTION(+:rho,tau,current,sdens,sodens) DO nst=1,nstloc - ps4=psi(:,:,:,:,nst) + ps4=psi(:,:,:,:,nst) CALL tstep(isospin(globalindex(nst)),mxpact/2,ps4) CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), & - ps4,rho,tau,current,sdens,sodens) + ps4,rho,tau,current,sdens,sodens) ENDDO !$OMP END PARALLEL DO IF(tmpi) CALL collect_densities @@ -256,7 +256,7 @@ SUBROUTINE dynamichf sodens=0.5D0*sodens sdens=0.5D0*sdens ! compute mean field and add external field - CALL skyrme + CALL skyrme IF(text_timedep) CALL extfld(time+dt/2.0D0) ! Step 3: full time step ! reset densities @@ -269,25 +269,25 @@ SUBROUTINE dynamichf !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,ps4) SCHEDULE(STATIC) & !$OMP REDUCTION(+:rho,tau,current,sdens,sodens) DO nst=1,nstloc - ps4=psi(:,:,:,:,nst) + ps4=psi(:,:,:,:,nst) CALL tstep(isospin(globalindex(nst)),mxpact,ps4) CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), & - ps4,rho,tau,current,sdens,sodens) + ps4,rho,tau,current,sdens,sodens) psi(:,:,:,:,nst)=ps4 ENDDO !$OMP END PARALLEL DO ! sum up over nodes IF(tmpi) CALL collect_densities ! Step 4: eliminate center-of-mass motion if desired - IF(mrescm/=0) THEN - IF(MOD(iter,mrescm)==0) THEN + IF(mrescm/=0) THEN + IF(MOD(iter,mrescm)==0) THEN CALL resetcm !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst) SCHEDULE(STATIC) & !$OMP REDUCTION(+:rho,tau,current,sdens,sodens) DO nst=1,nstloc CALL add_density(isospin(globalindex(nst)), & wocc(globalindex(nst)), & - psi(:,:,:,:,nst),rho,tau,current,sdens,sodens) + psi(:,:,:,:,nst),rho,tau,current,sdens,sodens) ENDDO !$OMP END PARALLEL DO IF(tmpi) CALL collect_densities @@ -298,24 +298,24 @@ SUBROUTINE dynamichf CALL tinfo ! Step 6: finishing up ! compute densities, currents, potentials etc. * - CALL skyrme + CALL skyrme IF(text_timedep) CALL extfld(time+dt) - IF(MOD(iter,mrest)==0) THEN + IF(MOD(iter,mrest)==0) THEN CALL write_wavefunctions IF(wflag) WRITE(*,*) ' Wrote restart file at end of iter=',iter ENDIF END DO Timestepping DEALLOCATE(ps4) END SUBROUTINE dynamichf -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: tstep !> @brief !!In this subroutine one wave function given as the argument \c psout -!!is stepped forward in time by the interval \c dt. +!!is stepped forward in time by the interval \c dt. !> !> @details !!The method used is the expansion of the exponential time-development operator -!!cut off at the power of \c mxp, which in practice is usually around 6. +!!cut off at the power of \c mxp, which in practice is usually around 6. !!Suppressing the argument of \f$ \hat h \f$ for brevity, we can write !!\f[ \hat U(t,t+\Delta t)\,\phi\approx \sum_{n=0}^m \phi^{(n)} \f] !!with @@ -335,25 +335,25 @@ END SUBROUTINE dynamichf !> REAL(db), takes the power to which the time-development operator is expanded. !> @param[in,out] psout !> REAL(db), array, takes the wave function and returns evolved wave function. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE tstep(iq,mxp,psout) INTEGER,INTENT(IN) :: iq,mxp COMPLEX(db),INTENT(INOUT) :: psout(:,:,:,:) INTEGER :: m,is,ix,iy,iz - COMPLEX(db) :: ps1(nx,ny,nz,2),ps2(nx,ny,nz,2) + COMPLEX(db) :: ps1(nx,ny,nz,2),ps2(nx,ny,nz,2) REAL(db) :: fmd ps1=psout !*********************************************************************** ! compute exp(-I*dt*h) by power series expansion * !*********************************************************************** DO m=1,mxp - fmd=-dt/(hbc*m) + fmd=-dt/(hbc*m) CALL hpsi(iq,esf,ps1,ps2) - DO is=1,2 - DO iz=1,nz - DO iy=1,ny - DO ix=1,nx - ps1(ix,iy,iz,is)=& + DO is=1,2 + DO iz=1,nz + DO iy=1,ny + DO ix=1,nx + ps1(ix,iy,iz,is)=& CMPLX(-fmd*AIMAG(ps2(ix,iy,iz,is)), & fmd*REAL(ps2(ix,iy,iz,is)),db) psout(ix,iy,iz,is)=psout(ix,iy,iz,is)+ps1(ix,iy,iz,is) @@ -363,11 +363,11 @@ SUBROUTINE tstep(iq,mxp,psout) ENDDO END DO END SUBROUTINE tstep - !--------------------------------------------------------------------------- + !--------------------------------------------------------------------------- ! DESCRIPTION: tinfo !> @brief !!This subroutine is used to output various pieces of information -!!relevant especially to the dynamic mode of the code. +!!relevant especially to the dynamic mode of the code. !> !> @details !!It is called at @@ -378,7 +378,7 @@ END SUBROUTINE tstep !! - Step 2: twobody analysis the twobody analysis is !! performed, but only if the calculation started as a twobody !! scenario. This is the simplified version which essentially only -!! determines the fragment in case of separation. +!! determines the fragment in case of separation. !! - Step 3: moments: the moments of the distribution are !! calculated using subroutine \c moments. This includes total mass, !! momenta, and angular momenta. They are printed out if indicated by @@ -392,9 +392,9 @@ END SUBROUTINE tstep !! because some of the logic may depend on them, especially the !! twobody-analysis, which needs the correct c.m., for example and is !! calculated at every time step and on every node. -!! +!! !! - Step 4: single-particle quantities: the single-particle -!! energies are calculated straightforwardly as expectation values +!! energies are calculated straightforwardly as expectation values !! of the Hamiltonian. !! The routine \c sp_properties is then called to obtain the other !! single-particle properties like angular momenta. They are @@ -411,13 +411,13 @@ END SUBROUTINE tstep !! The energies are protocolled in \c energiesfile and on standard !! output. !! - Step 6: density output: at intervals of \c mprint or in -!! the first time step the density printer plot is generated using +!! the first time step the density printer plot is generated using !! \c plot_densities and the binary densities are written onto *.tdd !! files using \c write_densities. !! - Step 7: other output: in the proper \c mprint interval !! the single-particle state !! information, and the moments are printed on standard output, using -!! also the routine \c moment_print. +!! also the routine \c moment_print. !! !! - Step 8: check for final separation: for the twobody case !! it is checked whether the separation found between the two fragments @@ -425,13 +425,13 @@ END SUBROUTINE tstep !! derivative \c rdot of the separation distance, in which case the !! program complete twobody analysis is done and the wave functions !! are saved; then the job is terminated with an appropriate message. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE tinfo REAL(db), DIMENSION(2) :: ecoll ! storage for collective-flow energy INTEGER :: il CHARACTER(*),PARAMETER :: & header=' # v**2 Norm Ekin Energy & - & Lx Ly Lz Sx Sy Sz ' + & Lx Ly Lz Sx Sy Sz ' !*********************************************************************** ! calculates dynamic observables for printout * !*********************************************************************** @@ -442,26 +442,26 @@ SUBROUTINE tinfo ! Step 1 printnow=mprint>0.AND.MOD(iter,mprint)==0 ! Step 2: twobody analysis - IF(nof/=2) THEN + IF(nof/=2) THEN istwobody=.FALSE. ELSE CALL twobody_analysis(.FALSE.) ! get distance in separated case ENDIF ! Step 3: moments calculated - + CALL moments(L_val,M_val) IF(printnow.AND.wflag) THEN - ! OPEN(unit=scratch,file=dipolesfile,POSITION='APPEND') + ! OPEN(unit=scratch,file=dipolesfile,POSITION='APPEND') ! WRITE(scratch,'(1x,i5,6E14.4)') iter,cmtot,cm(:,2)-cm(:,1) ! CLOSE(unit=scratch) - OPEN(unit=scratch,file=momentafile, POSITION='APPEND') + OPEN(unit=scratch,file=momentafile, POSITION='APPEND') WRITE(scratch,'(1x,f10.2,3g14.6)') time,pcm(:,1)+pcm(:,2) CLOSE(unit=scratch) CALL moment_shortprint() IF(texternal) CALL print_extfield() ENDIF ! Step 4: single-particle properties - IF(printnow) THEN + IF(printnow) THEN sp_energy=0.0D0 sp_norm=0.0D0 DO nst=1,nstloc @@ -485,10 +485,10 @@ SUBROUTINE tinfo (current(:,:,:,1,iq)**2+current(:,:,:,2,iq)**2& +current(:,:,:,3,iq)**2)/rho(:,:,:,iq) ) ENDDO - OPEN(unit=scratch,file=spinfile, POSITION='APPEND') - WRITE(scratch,'(1x,i5,9F10.4)') iter,orbital,spin,total_angmom + OPEN(unit=scratch,file=spinfile, POSITION='APPEND') + WRITE(scratch,'(1x,i5,9F10.4)') iter,orbital,spin,total_angmom CLOSE(unit=scratch) - OPEN(unit=scratch,file=energiesfile,POSITION='APPEND') + OPEN(unit=scratch,file=energiesfile,POSITION='APPEND') WRITE(scratch,'(F10.2,2F8.3,2F15.7,3(1PG13.5))') & time,pnr,ehf,ehfint,tke,ecoll CLOSE(unit=scratch) @@ -514,7 +514,7 @@ SUBROUTINE tinfo WRITE(*,'(/A)') ' Neutron Single Particle States:',header DO il=1,nstmax IF(il==npmin(2)) THEN - WRITE(*,'(/A)') ' Proton Single Particle States:',header + WRITE(*,'(/A)') ' Proton Single Particle States:',header END IF WRITE(*,'(1X,I3,F8.5,F9.6,F8.3,F10.3,3F8.3,3F7.3)') & il,wocc(il),sp_norm(il),sp_kinetic(il), & @@ -525,19 +525,19 @@ SUBROUTINE tinfo DEALLOCATE(ps1) ! Step 8: check whether final distance is reached and it is ! increasing in the twobody case - IF(istwobody.AND.roft>rsep.AND.roft>roft_old) THEN + IF(istwobody.AND.roft>rsep.AND.roft>roft_old) THEN CALL twobody_analysis(.TRUE.) ! do complete version CALL write_wavefunctions CALL write_densities IF(wflag) WRITE(*,*) ' Final separation distance reached' - STOP ' Final distance reached' + STOP ' Final distance reached' ENDIF initialcall=.FALSE. END SUBROUTINE tinfo -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: resetcm !> @brief -!!This subroutine resets the center-of-mass velocity to zero. +!!This subroutine resets the center-of-mass velocity to zero. !> !> @details !!The velocity, or rather the corresponding wave vector, is calculated from @@ -549,7 +549,7 @@ END SUBROUTINE tinfo !!\f] !!The wave functions are then multiplied by a common plane-wave phase !!factor \f$ \exp(-\I\vec k\cdot \vec r) \f$ to give a counter boost. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE resetcm INTEGER :: ix,iy,iz,is,nst REAL(db) :: akf(3),totmass @@ -560,7 +560,7 @@ SUBROUTINE resetcm END DO ! now apply correction to each wavefunction DO nst=1,nstloc - FORALL(is=1:2,ix=1:nx,iy=1:ny,iz=1:nz) + FORALL(is=1:2,ix=1:nx,iy=1:ny,iz=1:nz) psi(ix,iy,iz,is,nst)=& psi(ix,iy,iz,is,nst) * EXP( & CMPLX(0.D0,akf(1)*x(ix)+akf(2)*y(iy)+akf(3)*z(iz),db)) diff --git a/Code/energies.f90 b/Code/energies.f90 index 40efc99..34e08be 100644 --- a/Code/energies.f90 +++ b/Code/energies.f90 @@ -1,17 +1,17 @@ !------------------------------------------------------------------------------ ! MODULE: Modulename !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module computes the total energy and the various contributions to -!!it in two ways. The first method evaluates the density functional, +!!it in two ways. The first method evaluates the density functional, !!by direct integration to compute the integrated !!energy \c ehfint. The second method uses the sum of !!single-particle energies plus the rearrangement energies, \f$ E_{3,\rm !!corr} \f$ for the density dependent part and \f$ E_{C,\rm corr} \f$ for !!Coulomb exchange. !> -!>@details +!>@details !!The two ways of calculating the energy are assigned to the subroutines !!\c integ_energy, which also calculates the rearrangement energies, !!and \c sum_energy. Note that since \c integ_energy is always @@ -53,7 +53,7 @@ !! contact pairing interaction involving the pairing densities \f$ \xi_q \f$ !! augmented by an optional density dependence. The formula is !! \f[ E_{\rm pair} = \frac{1}{4} \sum_{q\in\{p,n\}}V_\mathrm{pair,q} -!! \int \D^3r |\xi_q|^2 \left[1 -\frac{\rho}{\rho_{0,\mathrm{pair}}}\right]\;. \f] +!! \int \D^3r |\xi_q|^2 \left[1 -\frac{\rho}{\rho_{0,\mathrm{pair}}}\right]\;. \f] !! It contains a continuous switch, the parameter !! \f$ \rho_{0,\mathrm{pair}} \f$. A pure \f$ \delta \f$-interaction (DI), also !! called volume pairing, is recovered for @@ -101,7 +101,7 @@ MODULE Energies REAL(db) :: total_angmom(3) !< the three components of the total !!angular momentum in units of \f$ \hbar \f$. CONTAINS -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: integ_energy !> @brief !!The purpose of this subroutine is to calculate the integrated energy. @@ -118,28 +118,28 @@ MODULE Energies !! !!The calculation proceeds in the following steps: !! - Step 1: the Laplacian of the densities is calculated in -!! \c worka, then the integrals for +!! \c worka, then the integrals for !! \c ehf0,\c ehf2, and \c ehf3 are performed. !! After the loop the result for \c ehf3 is also used to calculate !! \c e3corr. !! - Step 2: the integral for \c ehf1 is evaluated !! using \c worka for the \f$ \vec\jmath_q{}^2 \f$ term. -!! - Step 3: the spin-orbit contribution of +!! - Step 3: the spin-orbit contribution of !! \c ehfls is calculated using \c worka as storage for !! \f$ \nabla\cdot\vec J_q \f$. -!! - Step 4: the Coulomb energy \c ehfc is evaluated +!! - Step 4: the Coulomb energy \c ehfc is evaluated !! with the Slater correction taken into account if the !! force's \c ex is nonzero. At the same time the Coulomb correction !! for the summed energy is calculated !! and stored in \c ecorc. It will be used in !! the subroutine \c sum_energy. -!! - Step 5: the kinetic energy is integrated for +!! - Step 5: the kinetic energy is integrated for !! \c ehft. Note that only at this point the correct prefactor !! \f$ \hbar^2/2m \f$ is added; the use of \c tau in other expressions !! assumes its absence. !! - Step 6: Finally all terms are added to produce the total !! energy, \c efundet, from which the pairing energies are subtracted. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE integ_energy USE Trivial, ONLY: rmulx,rmuly,rmulz USE Grids, ONLY: wxyz,der1x,der2x,der1y,der2y,der1z,der2z @@ -234,27 +234,27 @@ SUBROUTINE integ_energy ! Step 6: form total energy ehfint=ehft+ehf0+ehf1+ehf2+ehf3+ehfls+ehfc-epair(1)-epair(2) END SUBROUTINE integ_energy -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: sum_energy !> @brief !!This subroutine mainly computes the Koopman sum, but also -!!sums up a number of other single-particle properties. +!!sums up a number of other single-particle properties. !> !> @details !!For systematics, the latter should be done in a different !!place, but at present is left here. !! -!!The summation of the total energy uses \c spenerg to compute +!!The summation of the total energy uses \c spenerg to compute !!\f[ \sum_k (\epsilon_k-\tfrac1{2}v_k)=\tfrac1{2}\sum_k(2t_k+v_k)= !!\tfrac1{2}\sum_k(t_k+\epsilon_k). \f] -!!The last sum is calculated, the rearrangement corrections +!!The last sum is calculated, the rearrangement corrections !!are added and the pairing energies subtracted. !! !!The subroutine then sums up the single-particle energy fluctuation !!\c sp_efluct1 and \c sp_efluct2, dividing them by the nucleon !!number. Finally the orbital and spin angular -!!momentum components are summed to form the total ones. -!--------------------------------------------------------------------------- +!!momentum components are summed to form the total ones. +!--------------------------------------------------------------------------- SUBROUTINE sum_energy USE Moment, ONLY: pnrtot INTEGER :: i diff --git a/Code/external.f90 b/Code/external.f90 index 1214bd1..17a1f9e 100644 --- a/Code/external.f90 +++ b/Code/external.f90 @@ -1,21 +1,21 @@ !------------------------------------------------------------------------------ ! MODULE: External !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module allows the coupling of the nucleonic wave functions to an -!!external field. +!!external field. !> -!>@details +!>@details !!This can be done either by adding a time-dependent external (i.e., not -!!self-consistent) potential to the single-particle Hamiltonian. +!!self-consistent) potential to the single-particle Hamiltonian. !!The time-dependence can be of Gaussian form !!\f[ f(t)=\exp\left(-(t-\tau_0)^2/\Delta\tau^2\right)\cos(\omega(\tau-\tau_0)) \f] !!or cosine squared form !!\f[ f(t)=\cos\left(\frac\pi{2}\left(\frac{t-\tau_0}{\Delta\tau}\right)^2\right) !! \theta\left(\Delta\tau-|t-\tau_0|\right)\cos(\omega(\tau-\tau_0)) \f] !!or by -!!giving an initial "boost" +!!giving an initial "boost" !!\f[ \psi_k(\vec r,s,t\!=\!0)= !!\psi_{k,0}(\vec r,s)\,\exp\left(-\I \eta F_q(\vec r)\right) \f] !!to each wave function. Since this is very @@ -25,16 +25,16 @@ !!assumptions and the isospin are however, fully functional and should !!be useful in many cases. !! -!!The operator \f$ F_q \f$ is the excitation operator which can be of any form from monopole +!!The operator \f$ F_q \f$ is the excitation operator which can be of any form from monopole !! to dipole, quadrupole, octupole and so on. However to avoid the spurious states in the strength !! the definition in the case of isoscalar dipole (L=1) is different than all the other as !!\f[ F_{1M} = (r^3 - \frac{5}{3} r)Y_{1M} \f] -!! where \f$ \f$ is the average value of \f$r^2\f$ which is equal to \f$r_{rms}^2\f$ and +!! where \f$ \f$ is the average value of \f$r^2\f$ which is equal to \f$r_{rms}^2\f$ and !!\f$r_{rms}\f$ is given as input variable r_avg for this case. !! For all other cases \f$ F_q \f$ is defined as !! \f[ F_{LM} = \sqrt{2L+1} r^L Y_{LM} \f] !! where \f$ Y_{LM} \f$ are the spherical harmonics. Only damped version of boundary conditions is implemented -!! in the present work. +!! in the present work. !------------------------------------------------------------------------------ MODULE External USE Params, ONLY: db,pi,wflag,extfieldfile,time,scratch @@ -43,9 +43,9 @@ MODULE External USE Levels, ONLY: nstloc,isospin,charge_number,mass_number USE MEANFIELD, ONLY: upot Use Spherical_Harmonics - IMPLICIT NONE + IMPLICIT NONE INTEGER :: isoext=0 !< isospin behavior of the external field: \c isoext - !! denotes the same action on protons and neutrons, + !! denotes the same action on protons and neutrons, !! isoext=1 that with opposing signs. INTEGER,PRIVATE :: ipulse=0 !< type of pulse. ipulse=0 denotes the !! initial boost configuration, ipulse=1 a Gaussian pulse @@ -61,43 +61,43 @@ MODULE External REAL(db),PRIVATE :: taut !< parameter \f$ \Delta\tau \f$ for time-dependent excitation REAL(db),PRIVATE :: tau0 !< parameter \f$ \tau_0 \f$ for time-dependent excitation LOGICAL,PRIVATE :: textfield_periodic=.false. !< if this is set to \c true, the - !! external field is made periodic by substituting - !! \f$ x^2\rightarrow\sin^2\left(\pi x/x_L\right) \f$, + !! external field is made periodic by substituting + !! \f$ x^2\rightarrow\sin^2\left(\pi x/x_L\right) \f$, !! otherwise a damping factor is used \f$ F_q(\vec r) !! \rightarrow \frac{F_q(\vec r)}{1+e^{(r-r_0)/\Delta r}} \f$. REAL(db) :: ampl_ext=0.D0 !< The strength parameter (\f$ \eta \f$) for the perturbation. - !! Its numerical magnitude is usually not important by itself, - !! but varying it allows studying the effects of different + !! Its numerical magnitude is usually not important by itself, + !! but varying it allows studying the effects of different !! strengths of the excitation. INTEGER :: L_val=0 !< paramter \f$L_{val}\f$, the orbital angular quantum number(L) used - !! in the definition of spherical harmoninc (\f$Y_{L}^{M}\f$) + !! in the definition of spherical harmoninc (\f$Y_{L}^{M}\f$) !! which is used while defining the external field INTEGER :: M_val=0 !< paramter \f$M_{val}\f$, the projection quantum number(M) used - !! in the definition of spherical harmoninc (\f$Y_{L}^{M}\f$) + !! in the definition of spherical harmoninc (\f$Y_{L}^{M}\f$) !! which is used while defining the external field - INTEGER,PRIVATE :: only_P=0 !< If this is set to 1, then the strength is calculated by + INTEGER,PRIVATE :: only_P=0 !< If this is set to 1, then the strength is calculated by !! only exciting the protons to calculate EM response. It should be used carefully - !! because there is no COM correction implemented for this, however, it works fine in the - !! cases of big doubly magic nuclei like \f$^{208}Pb\f$. + !! because there is no COM correction implemented for this, however, it works fine in the + !! cases of big doubly magic nuclei like \f$^{208}Pb\f$. REAL(db) :: r_avg=0.0d0 !< The value of \f$ \sqrt{} \f$ used in the case of dipole boost. PUBLIC ::L_val,M_val,ampl_ext,isoext SAVE CONTAINS -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: getin_external !> @brief !!It reads in all the parameters of the external field. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE getin_external NAMELIST/extern/ ampl_ext,L_val,M_val,radext,widext,isoext,& ipulse,omega,tau0,taut,textfield_periodic,r_avg,only_P READ(5,extern) END SUBROUTINE getin_external -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: init_external !> @brief -!!It does some consistency checks of the external parameters and initializes the external field. +!!It does some consistency checks of the external parameters and initializes the external field. !> !> @details !!The relative strength for the neutron and proton fields is set equal for @@ -108,8 +108,8 @@ END SUBROUTINE getin_external !!Then the array \c extfield is allocated and the time-independent !!spatial potential calculated for both isospin cases and in either the !! Here only the damped version of the boundary conditions -!! are implemented for the external field \f$F_q\f$. -!--------------------------------------------------------------------------- +!! are implemented for the external field \f$F_q\f$. +!--------------------------------------------------------------------------- SUBROUTINE init_external REAL(db) :: facn,facp,facr,xlim,ylim,zlim,dip_f INTEGER :: ix,iy,iz @@ -118,8 +118,8 @@ SUBROUTINE init_external IF(ipulse<0.OR.ipulse>2) STOP & ' External field: called with invalid pulse type' IF(wflag) THEN - WRITE(*,*) "***** Parameters of external field *****" - WRITE(*,"(a,e12.4)") " Amplitude of axial quad. =",ampl_ext + WRITE(*,*) "***** Parameters of external field *****" + WRITE(*,"(a,e12.4)") " Amplitude of axial quad. =",ampl_ext WRITE(*,"(2(A,F10.4),A)") " Radial damping: radius ",radext, & ' fm, width ',widext,' fm' WRITE(*,"(2(A,I2),2A)") " Isospin of excitation:",isoext, & @@ -136,10 +136,10 @@ SUBROUTINE init_external ' External field: tau0 @brief !!This is again a very straightforward routine. It calculates the @@ -204,19 +204,19 @@ END SUBROUTINE init_external !> !> @param[in] time !> REAL(db), takes the current time. -!--------------------------------------------------------------------------- - SUBROUTINE extfld(time) +!--------------------------------------------------------------------------- + SUBROUTINE extfld(time) REAL(db),INTENT(IN) :: time REAL(db) :: time_factor - IF(ipulse==1) THEN - time_factor=EXP(-((time-tau0)/taut) **2) - ELSE - IF(timetau0+taut) THEN - time_factor=0.0D0 - ELSE - time_factor=COS(0.5D0*pi *(time-tau0)/taut) **2 + IF(ipulse==1) THEN + time_factor=EXP(-((time-tau0)/taut) **2) + ELSE + IF(timetau0+taut) THEN + time_factor=0.0D0 + ELSE + time_factor=COS(0.5D0*pi *(time-tau0)/taut) **2 ENDIF ENDIF IF(omega/=0.0D00) THEN @@ -224,10 +224,10 @@ SUBROUTINE extfld(time) ENDIF upot=time_factor*extfield + upot END SUBROUTINE extfld -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: extboost !> @brief -!!This performs the initial boost on all single-particle wave functions. +!!This performs the initial boost on all single-particle wave functions. !> !> @details !!Note that it is always called @@ -240,7 +240,7 @@ END SUBROUTINE extfld !> !> @param[out] noboostflag !> LOGICAL, returns if boost has to be performed. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE extboost(noboostflag) USE Levels, ONLY: psi LOGICAL,INTENT(OUT) :: noboostflag @@ -253,7 +253,7 @@ SUBROUTINE extboost(noboostflag) -extfield(:,:,:,isospin(globalindex(nst))),db)) END FORALL END SUBROUTINE extboost -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: print_extfield !> @brief !!This subroutine calculates the expectation of the coupling energy to @@ -261,10 +261,10 @@ END SUBROUTINE extboost !!\f[ \sum_q\int\,\D^3r\,\rho_q(\vec r)F_q(\vec r) \f] !!and prints one line containing the present time and this value onto !!the file \c extfieldfile. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE print_extfield() USE Densities, ONLY: rho - OPEN(UNIT=scratch,file=extfieldfile,POSITION='APPEND') + OPEN(UNIT=scratch,file=extfieldfile,POSITION='APPEND') WRITE(scratch,'(3x,F12.3,3x,F35.15,2x,F25.18,2I5)') time,wxyz*SUM(rho*extfield)/ampl_ext ! print*,'External',time,wxyz*SUM(rho*extfield)/ampl_ext CLOSE(UNIT=scratch) diff --git a/Code/forces.f90 b/Code/forces.f90 index 5f5d3fe..13b0983 100644 --- a/Code/forces.f90 +++ b/Code/forces.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! MODULE: Modulename !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!Module \c Forces describes the interactions used in the code. The !!idea is to produce a library of Skyrme forces that can be called up @@ -27,8 +27,8 @@ MODULE Forces REAL(db) :: v0prot! This contains the parameters for the Skyrme force: @@ -38,9 +38,9 @@ MODULE Forces !!exchange term. For ex=1 it is included (this is the normal !!case), for ex=0 not. INTEGER :: zpe !@name Skyrme parameters. + !>@name Skyrme parameters. !>@{ REAL(db) :: t0,t1,t2,t3,t4 !>@} @@ -49,7 +49,7 @@ MODULE Forces REAL(db) :: x0,x1,x2,x3,b4p !>@} REAL(db) :: power !< exponent in the nonlinear (originally three-body) term. - TYPE(Pairing) :: vdi !< parameter set for the volume-delta pairing case. + TYPE(Pairing) :: vdi !< parameter set for the volume-delta pairing case. TYPE(Pairing) :: dddi !< parameter set for the density-dependent delta pairing case. END TYPE Force ! include predefined forces @@ -59,7 +59,7 @@ MODULE Forces !!historical reasons the values are 0: no pairing, 5: VDI pairing, and !!6: DDDI pairing. In the input the symbolic names are used so these !!numerical values are hidden to the user. For details see the input - !!description and module \c Pairs. + !!description and module \c Pairs. TYPE(Force) :: f !< this contains parameters for !!the Skyrme force actually used in the present calculation, packed !into the derived-type \c Force. @@ -72,15 +72,15 @@ MODULE Forces REAL(db) :: nucleon_mass !< the mass of the nucleon (average of !!neutron and Proton) in MeV calculated from \c h2ma and \c hbc. ! - !>@name these are the coefficients actually used + !>@name these are the coefficients actually used !!for the mean-field and single-particle Hamiltonian calculations - !!in \c skyrme and \c integ_energy. Note that only \c b4p is also included + !!in \c skyrme and \c integ_energy. Note that only \c b4p is also included !!in the Skyrme-force definition; the others are derived from the \c t coefficients. !>@{ REAL(db) :: b0,b0p,b1,b1p,b2,b2p,b3,b3p,b4,b4p,slate !>@} CONTAINS -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: Routinename !> @brief !!The purpose of the subroutine is to read the force and pairing @@ -119,7 +119,7 @@ MODULE Forces !!Finally the routine calculates the values of \c nucleon_mass and !!\c h2ma. It then prints out a description of the force and pairing !!parameters. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE read_force CHARACTER(8) :: name,pairing INTEGER :: ex,zpe @@ -136,7 +136,7 @@ SUBROUTINE read_force ex,zpe,h2m,t0,t1,t2,t3,t4,x0,x1,x2,x3,b4p,power, & ipair,v0prot,v0neut,rho0pr,mixture,turnoff_zpe ! mark force & pairing parameters as undefined - h2m=-1.0; v0prot=-1.0; v0neut=-1.0; + h2m=-1.0; v0prot=-1.0; v0neut=-1.0; READ(5,force) ! seek for force in predefined ones predefined=.FALSE. @@ -168,16 +168,16 @@ SUBROUTINE read_force WRITE(*,*) '***** Zero-point-energy correction turned off' END IF ! calculate "b" and Slater coefficients - b0=f%t0*(1.0D0+0.5D0*f%x0) - b0p=f%t0*(0.5D0+f%x0) - b1=(f%t1+0.5D0*f%x1*f%t1+f%t2+0.5*f%x2*f%t2)/4.0D0 - b1p=(f%t1*(0.5D0+f%x1)-f%t2*(0.5D0+f%x2))/4.0D0 + b0=f%t0*(1.0D0+0.5D0*f%x0) + b0p=f%t0*(0.5D0+f%x0) + b1=(f%t1+0.5D0*f%x1*f%t1+f%t2+0.5*f%x2*f%t2)/4.0D0 + b1p=(f%t1*(0.5D0+f%x1)-f%t2*(0.5D0+f%x2))/4.0D0 b2=(3.0D0*f%t1*(1.D0+0.5D0*f%x1)-f%t2*(1.D0+0.5D0*f%x2))/8.0D0 b2p=(3.D0*f%t1*(0.5D0+f%x1)+f%t2*(0.5D0+f%x2))/8.D0 b3=f%t3*(1.D0+0.5D0*f%x3)/4.D0 - b3p=f%t3*(0.5D0+f%x3)/4.D0 - b4=f%t4/2.D0 - b4p=f%b4p + b3p=f%t3*(0.5D0+f%x3)/4.D0 + b4=f%t4/2.D0 + b4p=f%b4p slate=(3.0D0/pi)**(1.0D0/3.0D0)*e2 ! now set up pairing: first case of none IF(TRIM(pairing)=='NONE') THEN @@ -210,7 +210,7 @@ SUBROUTINE read_force WRITE(*,"(5(A6,F12.5))") "t0",f%t0,"t1",f%t1,"t2",f%t2,"t3",f%t3,"t4",f%t4 WRITE(*,"(5(A6,F12.5))") "x0",f%x0,"x1",f%x1,"x2",f%x2,"x3",f%x3,"b4p",f%b4p WRITE(*,"(A6,F12.5)") "Power",f%power - WRITE(*,"(A,I2)") " Pairing parameters: Option ipair:",ipair + WRITE(*,"(A,I2)") " Pairing parameters: Option ipair:",ipair WRITE(*,"(3(A7,F12.5))") "v0prot",p%v0prot,"v0neut",p%v0neut,"rho0pr",p%rho0pr,'MIX',p%mixture ENDIF END SUBROUTINE read_force diff --git a/Code/fourier.f90 b/Code/fourier.f90 index 950dc41..fcffb7c 100644 --- a/Code/fourier.f90 +++ b/Code/fourier.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! MODULE: Fourier !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module initializes the \c FFTW3 package for doing Fourier !!transforms \cite Fri05a. It needs the file \c fftw3.f from the @@ -41,7 +41,7 @@ MODULE Fourier INTEGER(C_LONG),SAVE :: zforward,zbackward !>@} CONTAINS -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: init_fft !> @brief !!In this subroutine the \c FFTW system is initialized @@ -67,11 +67,11 @@ MODULE Fourier !!these calls: in the code all transforms are in-place except for \c xforward, !!\c yforward, and \c zforward. It was found that !!very strange things happen if this rule is not obeyed. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE init_fft INCLUDE 'fftw3.f' COMPLEX(db),ALLOCATABLE :: p(:,:,:,:,:) - INTEGER,SAVE :: FFTW_planflag + INTEGER,SAVE :: FFTW_planflag ! set option for FFTW setup here ! FFTW_planflag=FFTW_ESTIMATE ! FFTW_planflag=FFTW_MEASURE diff --git a/Code/fragments.f90 b/Code/fragments.f90 index 2f338bb..b8b818a 100644 --- a/Code/fragments.f90 +++ b/Code/fragments.f90 @@ -1,13 +1,13 @@ !------------------------------------------------------------------------------ ! MODULE: Fragments !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module is concerned with setting up the initial condition from !!data files containing fragment wave functions, usually obtained in a -!!previous static calculation. +!!previous static calculation. !> -!>@details +!>@details !!A typical application is the initialization of a heavy-ion reaction. !!Beyond that, any number of fragments (limited by the parameter !!variable \c mnof) can be put into the grid at prescribed positions @@ -29,10 +29,10 @@ !------------------------------------------------------------------------------ MODULE Fragments USE Params - USE Grids + USE Grids USE Forces, ONLY: f,nucleon_mass USE Levels - IMPLICIT NONE + IMPLICIT NONE SAVE PRIVATE LOGICAL :: fix_boost ! @brief !!This subroutine reads the input variables, makes some consistency @@ -124,8 +124,8 @@ MODULE Fragments !! !!At the end the two-body initialization \c twobody_init is called !!for the dynamic case with two fragments and fix_boost=.FALSE.. -!!This calculates the \c fboost values from \c ecm and \c b. -!--------------------------------------------------------------------------- +!!This calculates the \c fboost values from \c ecm and \c b. +!--------------------------------------------------------------------------- SUBROUTINE getin_fragments INTEGER :: i REAL(db) :: fdx,fdy,fdz,fwxyz @@ -141,9 +141,9 @@ SUBROUTINE getin_fragments fix_boost=.FALSE. READ(5,fragments) END IF - IF(nof/=2.AND..NOT.fix_boost) THEN - IF(wflag) WRITE(*,"(//,a,//)") "Non-fixed boost only for nof=2" - STOP + IF(nof/=2.AND..NOT.fix_boost) THEN + IF(wflag) WRITE(*,"(//,a,//)") "Non-fixed boost only for nof=2" + STOP ENDIF IF(wflag) THEN IF(fix_boost) THEN @@ -177,7 +177,7 @@ SUBROUTINE getin_fragments filename(i) WRITE(*,"(A,1P,3E12.4)") " Location(x,y,z): ",fcent(:,i) WRITE(*,"(2(A,F9.4))") " Mass number: ",fmass(i), & - ' Charge number: ',fcharge(i) + ' Charge number: ',fcharge(i) IF(fix_boost) & WRITE(*,"(A,1P,3E12.4)") " Boost(x,y,z): ",fboost(:,i) END IF @@ -208,27 +208,27 @@ SUBROUTINE getin_fragments mass_number=nneut+nprot IF(.NOT.fix_boost.AND.nof==2.AND.tdynamic) CALL twobody_init END SUBROUTINE getin_fragments -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: read_fragments !> @brief !!This subroutine prints summary information about which fragment wave !!functions occupy which range of indices. Then it does a loop over -!!fragments to read in their wave functions using +!!fragments to read in their wave functions using !!\c read_one_fragment, followed by applying the boost to them using !!\c boost_fragment. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE read_fragments INTEGER :: iff IF(wflag) & WRITE(*,*) '***** Input of fragment wave functions *****' - DO iff=1,nof + DO iff=1,nof CALL read_one_fragment(iff) - CALL boost_fragment(iff) + CALL boost_fragment(iff) ENDDO IF(wflag) & WRITE(*,*) '***** All fragments loaded and boosted *****' END SUBROUTINE read_fragments -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: read_one_fragment !> @brief !!This subroutine reopens a fragment file @@ -245,13 +245,13 @@ END SUBROUTINE read_fragments !!this information is printed. !! !!Once this loop is concluded, the spatial shift is prepared. The shift -!!is calculated from the difference between the desired position \c fcent +!!is calculated from the difference between the desired position \c fcent !!with respect to the origin of the new coordinates given by !!\c x etc., and the fragment center-of-mass \c fcmtot with !!respect to the origin in fragment coordinates \c fx etc. !!Subroutine \c phases is used to calculate essentially the shift !!phase factor \f$ \exp(-\I\vec k\cdot\Delta\vec r) \f$, which is a product of -!!phases in each coordinate direction \c akx, \c aky, and \c akz. +!!phases in each coordinate direction \c akx, \c aky, and \c akz. !!These have an index corresponding to \f$ \vec k \f$ in the Fourier transform. !! !!Now the index arrays \c fnode and \c flocalindex are read, @@ -265,7 +265,7 @@ END SUBROUTINE read_fragments !!except for ignoring the input record. !! !!For the case of multiple files for one fragment (which is recognized -!!by not all the indices \c fnode being zero, a short subroutine \c locate +!!by not all the indices \c fnode being zero, a short subroutine \c locate !!is used to position input at the correct location. !! !!Otherwise the wave function is read from the fragment file using the @@ -273,11 +273,11 @@ END SUBROUTINE read_fragments !!full new dimension, then it is Fourier transformed, multiplied by the !!phase factor, and transformed back. Finally it is inserted into the !!total wave function array \c psi, where any zeroes are replaced by -!!a small number, presently set to \f$ 10^{-20} \f$. +!!a small number, presently set to \f$ 10^{-20} \f$. !> !> @param[in] iff !> INTEGER, takes the umber of the fragment which is read in. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE read_one_fragment(iff) USE Parallel, ONLY: node,mpi_myproc,localindex USE Fourier @@ -293,8 +293,8 @@ SUBROUTINE read_one_fragment(iff) INTEGER,DIMENSION(fnstmax(iff)) :: fnode,flocalindex OPEN(UNIT=scratch,FILE=filename(iff),STATUS='old',FORM=& 'unformatted') - READ(scratch) - READ(scratch) + READ(scratch) + READ(scratch) READ(scratch) fx,fy,fz ! read the s.p.quantities, reinserting them into the correct positions READ(scratch) fwocc,fsp_energy,fsp_parity,fsp_norm,fsp_kinetic,fsp_efluct1 @@ -302,7 +302,7 @@ SUBROUTINE read_one_fragment(iff) IF(wflag) THEN WRITE(*,'(A,I3)') ' ***** S.p. levels for fragment #',iff WRITE(*,'(A,3I5,A,2I5)') ' Iso #new #old Occup. E_sp & - &Parity Norm E_kin E_fluct' + &Parity Norm E_kin E_fluct' END IF DO iq=1,2 DO inew=fnewnpmin(iq,iff),fnewnpsi(iq,iff) @@ -344,7 +344,7 @@ SUBROUTINE read_one_fragment(iff) DO nst=1,fnumber(iq,iff) newnst=fnewnpmin(iq,iff)+nst-1 oldnst=fnpmin(iq,iff)+nst-1 ! includes empty states - IF(node(newnst)==mpi_myproc) THEN + IF(node(newnst)==mpi_myproc) THEN ipn=localindex(newnst) ! open correct file and position for multifile case IF(multifile) CALL locate(fnode(oldnst),flocalindex(oldnst)) @@ -352,12 +352,12 @@ SUBROUTINE read_one_fragment(iff) READ(scratch) ps1(1:fnx(iff),1:fny(iff),1:fnz(iff),:) DO is=1,2 CALL dfftw_execute_dft(pforward,ps1(:,:,:,is),ps1(:,:,:,is)) - FORALL(ix=1:nx,iy=1:ny,iz=1:nz) + FORALL(ix=1:nx,iy=1:ny,iz=1:nz) ps1(ix,iy,iz,is)=ps1(ix,iy,iz,is)*akx(ix)*aky(iy)*akz(iz) & /DBLE(nx*ny*nz) END FORALL CALL dfftw_execute_dft(pbackward,ps1(:,:,:,is),ps1(:,:,:,is)) - WHERE(ABS(ps1)==0.D0) + WHERE(ABS(ps1)==0.D0) psi(:,:,:,:,ipn)=1.0D-20 ELSEWHERE psi(:,:,:,:,ipn)=ps1 @@ -370,7 +370,7 @@ SUBROUTINE read_one_fragment(iff) ENDDO CLOSE(unit=scratch) END SUBROUTINE read_one_fragment -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: Routinename !> @brief !!This has the task to position the file for reading wave functions at @@ -394,7 +394,7 @@ END SUBROUTINE read_one_fragment !> INTEGER, takes the number of the file. !> @param[in] pos !> INTEGER, takes the index of the wave function in this file. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! position scratch file in correct location SUBROUTINE locate(fileno,pos) INTEGER,INTENT(IN) :: fileno,pos @@ -410,7 +410,7 @@ SUBROUTINE locate(fileno,pos) END DO END IF END SUBROUTINE locate -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: Routinename !> @brief !!This subroutine calculates the phase factors for a one-dimensional @@ -428,8 +428,8 @@ END SUBROUTINE locate !> @param[out] a !> COMPLEX(db), array, returns the phases. !> @param[in] c -!> REAL(db), takes the shift -!--------------------------------------------------------------------------- +!> REAL(db), takes the shift +!--------------------------------------------------------------------------- PURE SUBROUTINE phases(n,a,c) INTEGER :: n,i,k COMPLEX(db),INTENT(OUT) :: a(n) @@ -441,7 +441,7 @@ PURE SUBROUTINE phases(n,a,c) a(i)=EXP(CMPLX(0.D0,-2.D0*pi*k*c,db)) END DO END SUBROUTINE phases -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: twobody_init !> @brief !!The purpose of this subroutine is to calculate the boost values from @@ -458,7 +458,7 @@ END SUBROUTINE phases !! !!The Coulomb energy is calculated assuming two point charges at !!distance \c roft and is subtracted from \c ecm to yield the -!!kinetic energy remaining at this distance, from which the relative +!!kinetic energy remaining at this distance, from which the relative !!velocity \c vrel_d is calculated. Since the total center !!of mass is assumed to be at rest, the velocities of the fragments \c v1 !!and \c v2 can then be simply obtained. The instantaneous impact parameter @@ -471,8 +471,8 @@ END SUBROUTINE phases !!of \c fboost as given in the input. Note that these energies are !!signed to indicate the direction of motion. !! -!!Both velocities and energies are printed. -!--------------------------------------------------------------------------- +!!Both velocities and energies are printed. +!--------------------------------------------------------------------------- SUBROUTINE twobody_init REAL(db) :: vrel,vrel_d,b_d,v1,v2,ec,dix,diz,totmass,sint,cost, & xli,xmu,roft @@ -504,7 +504,7 @@ SUBROUTINE twobody_init WRITE(*,*) '***** Two-body initialization *****' WRITE(*,"(2(A,F12.4),A)") " c. m. Energy:",ecm, & ' MeV. Impact parameter b: ',b,' fm' - WRITE(*,"(A,F12.4)") " xli/hbar :",xli + WRITE(*,"(A,F12.4)") " xli/hbar :",xli WRITE(*,*) 'Computed boost velocities in units of c' WRITE(*,'(A,I2,3G15.6)') (' Fragment #',i,fboost(:,i),i=1,2) END IF @@ -518,12 +518,12 @@ SUBROUTINE twobody_init WRITE(*,'(A,I2,3G15.6)') (' Fragment #',i,fboost(:,i),i=1,2) END IF END SUBROUTINE twobody_init -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: boost_fragment !> @brief !!This subroutine multiplies the configuration-space wave functions of !!fragment no. \c iff by plane-wave factors to give them -!!translational motion. +!!translational motion. !> !> @details !!This is done by calculating the wave number @@ -539,23 +539,23 @@ END SUBROUTINE twobody_init !!The rest is then straightforward. !> !> @param[in] iff -!> INTEGER, takes the number of the fragment to be boosted. -!--------------------------------------------------------------------------- - SUBROUTINE boost_fragment(iff) +!> INTEGER, takes the number of the fragment to be boosted. +!--------------------------------------------------------------------------- + SUBROUTINE boost_fragment(iff) USE Parallel, ONLY: node,localindex,mpi_myproc INTEGER,INTENT(IN) :: iff REAL(db) :: sb(3),tb(3),akf(3) INTEGER :: iq,nst,is,ix,iy,iz - sb=0.0D0 - tb=fboost(:,iff) - WHERE(tb/=0.0D0) - sb=tb/ABS(tb) + sb=0.0D0 + tb=fboost(:,iff) + WHERE(tb/=0.0D0) + sb=tb/ABS(tb) END WHERE - DO iq=1,2 - akf=sb*SQRT(1.0D0/f%h2m(iq)*ABS(tb)/fmass(iff)) + DO iq=1,2 + akf=sb*SQRT(1.0D0/f%h2m(iq)*ABS(tb)/fmass(iff)) DO nst=fnewnpmin(iq,iff),fnewnpsi(iq,iff) IF(node(nst)==mpi_myproc) THEN - FORALL(is=1:2,ix=1:nx,iy=1:ny,iz=1:nz) + FORALL(is=1:2,ix=1:nx,iy=1:ny,iz=1:nz) psi(ix,iy,iz,is,localindex(nst))=& psi(ix,iy,iz,is,localindex(nst)) * EXP( & CMPLX(0.D0,akf(1)*x(ix)+akf(2)*y(iy)+akf(3)*z(iz),db)) @@ -564,7 +564,7 @@ SUBROUTINE boost_fragment(iff) ENDDO ENDDO END SUBROUTINE boost_fragment -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: errormsg !> @brief !!Writes out error messages. @@ -573,7 +573,7 @@ END SUBROUTINE boost_fragment !> CHARACTER, array, takes the first message. !> @param[in] msg2 !> CHARACTER, array, takes the second message. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE errormsg(msg1,msg2) CHARACTER(*),INTENT(IN) :: msg1,msg2 IF(wflag) WRITE(*,*) msg1,msg2 diff --git a/Code/grids.f90 b/Code/grids.f90 index c5d3826..db3c589 100644 --- a/Code/grids.f90 +++ b/Code/grids.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! MODULE: Modulename !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module deals with the definition of the spatial grid and associated operations. !------------------------------------------------------------------------------ @@ -15,43 +15,43 @@ MODULE Grids !!reflection symmetry. INTEGER :: nz !< Number of points in z-direction. Must be even to preserve !!reflection symmetry. - LOGICAL :: periodic !< logical variable indicating whether the situation is triply + LOGICAL :: periodic !< logical variable indicating whether the situation is triply !!periodic in three-dimensional space. REAL(db) :: dx !< The spacing between grid points (in fm) in x-directions. REAL(db) :: dy !< The spacing between grid points (in fm) in y-directions. REAL(db) :: dz !< The spacing between grid points (in fm) in z-directions. REAL(db) :: wxyz !< the volume element wxyz=dx*dy*dz. REAL(db),POINTER :: x(:) ! (nx,nx) and calculated in subroutine \c sder. - der2x, & ! (nx,nx) and calculated in subroutine \c sder2. - cdmpx, & ! (nx,nx) and calculated in subroutine \c setdmc. - der1y, & ! (ny,ny) and calculated in subroutine \c sder. - der2y, & ! (ny,ny) and calculated in subroutine \c sder2. - cdmpy, & ! (ny,ny) and calculated in subroutine \c setdmc. - der1z, & ! (nz,nz) and calculated in subroutine \c sder. - der2z, & ! (nz,nz) and calculated in subroutine \c sder2. - cdmpz ! (nz,nz) and calculated in subroutine \c setdmc. PRIVATE :: init_coord, sder, sder2, setdmc, gauss CONTAINS -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: init_grid !> @brief !!This subroutine is called during the initialization for both static @@ -63,8 +63,8 @@ MODULE Grids !> @details !!This is done by calling the subroutine \c init_coord once for each !!direction. It does everything needed except the calculation of -!!the volume element. -!--------------------------------------------------------------------------- +!!the volume element. +!--------------------------------------------------------------------------- SUBROUTINE init_grid NAMELIST /Grid/ nx,ny,nz,dx,dy,dz,periodic dx=0.D0 @@ -93,11 +93,11 @@ SUBROUTINE init_grid CALL init_coord('z',nz,dz,z,der1z,der2z,cdmpz) wxyz=dx*dy*dz END SUBROUTINE init_grid -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: Routinename !> @brief !!In this subroutine the defining information for a grid direction -!!(generically called \c v, which can be replaced by \c x, \c y, or \c z) +!!(generically called \c v, which can be replaced by \c x, \c y, or \c z) !!in the form of the number of points \c nv and the !!spacing \c dv is used to generate the associated data. The arrays !!of coordinate values, derivative and damping matrices are allocated. @@ -118,7 +118,7 @@ END SUBROUTINE init_grid !!be done by changing the statement generating the values of \c v. !! !!Finally the derivative matrices and damping matrix are computed using -!!\c sder, \c sder2, and \c setdmc. +!!\c sder, \c sder2, and \c setdmc. !> !> @param[in] name !> CHARACTER, array, takes x, y, or z as sirection. @@ -134,7 +134,7 @@ END SUBROUTINE init_grid !> REAL(db), array, returns matrix for second derivative. !> @param[out] cdmpv !> REAL(db), array, returns matrix for damping. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE init_coord(name,nv,dv,v,der1v,der2v,cdmpv) CHARACTER(*) :: name INTEGER :: nv @@ -152,17 +152,17 @@ SUBROUTINE init_coord(name,nv,dv,v,der1v,der2v,cdmpv) CALL sder(der1v,nv,dv) CALL sder2(der2v,nv,dv) END SUBROUTINE init_coord -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: sder !> @brief -!!This subroutine calculates the matrix for the first derivative. +!!This subroutine calculates the matrix for the first derivative. !> @param[out] der !> REAL(db), array, returns derivative matrix. !> @param[in] nmax !> INTEGER, takes the dimension. !> @param[in] d -!> REAL(db), takes the grid spacing. -!--------------------------------------------------------------------------- +!> REAL(db), takes the grid spacing. +!--------------------------------------------------------------------------- PURE SUBROUTINE sder(der,nmax,d) INTEGER :: nmax REAL(db) :: d,der(:,:) @@ -183,17 +183,17 @@ PURE SUBROUTINE sder(der,nmax,d) ENDDO ENDDO END SUBROUTINE sder -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: sder2 !> @brief -!!This subroutine calculates the matrix for the second derivative. +!!This subroutine calculates the matrix for the second derivative. !> @param[out] der !> REAL(db), array, returns derivative matrix. !> @param[in] nmax !> INTEGER, takes the dimension. !> @param[in] d -!> REAL(db), takes the grid spacing. -!--------------------------------------------------------------------------- +!> REAL(db), takes the grid spacing. +!--------------------------------------------------------------------------- PURE SUBROUTINE sder2(der,nmax,d) INTEGER :: nmax REAL(db) :: d,der(1:nmax,1:nmax) @@ -214,35 +214,35 @@ PURE SUBROUTINE sder2(der,nmax,d) ENDDO ENDDO END SUBROUTINE sder2 -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: setup_damping !> @brief !!This sets up the damping matrices by calls to \c setdmc for each -!!coordinate direction. +!!coordinate direction. !> !> @details -!!The reason for not including this in \c init_grid +!!The reason for not including this in \c init_grid !!is that it used only in the static calculation and !!requires the damping parameter \c e0dmp, which is in the static -!!module. It has to be passed as a parameter because +!!module. It has to be passed as a parameter because !!circular dependence of modules would result otherwise. !> !> @param[in] e0dmp !>REAL(db), takes the damping parameter. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE setup_damping(e0dmp) REAL(db),INTENT(IN) :: e0dmp CALL setdmc(der2x,nx,cdmpx,e0dmp) CALL setdmc(der2y,ny,cdmpy,e0dmp) CALL setdmc(der2z,nz,cdmpz,e0dmp) END SUBROUTINE setup_damping -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: setdmc !> @brief !!This subroutine calculates the matrices corresponding to the !!one-dimensional operators !!\f[ \frac{1}{1+\hat t/E_0} {\rm~with~} \hat -!!t=-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}. \f] +!!t=-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}. \f] !> !> @details !!Here \f$ E_0 \f$ is @@ -259,10 +259,10 @@ END SUBROUTINE setup_damping !> @param[in] nmax !> INTEGER, takes the number of grid points !> @param[out] cdmp -!> REAL(db), array, returns the damping matrix. +!> REAL(db), array, returns the damping matrix. !> @param[in] e0dmp !> REAL(db), takes the damping parameter. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- PURE SUBROUTINE setdmc(der2,nmax,cdmp,e0dmp) USE Forces, ONLY: h2ma REAL(db),INTENT(IN) :: der2(:,:) @@ -273,14 +273,14 @@ PURE SUBROUTINE setdmc(der2,nmax,cdmp,e0dmp) INTEGER :: i cdmp=0.D0 IF(e0dmp<=0.D0) RETURN - unit=-h2ma*der2/e0dmp + unit=-h2ma*der2/e0dmp FORALL(i=1:nmax) - unit(i,i)=1.0D0+unit(i,i) - cdmp(i,i)=1.0D0 + unit(i,i)=1.0D0+unit(i,i) + cdmp(i,i)=1.0D0 END FORALL CALL gauss(unit,cdmp,nmax) END SUBROUTINE setdmc -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: gauss !> @brief !!This is a Fortran 95 implementation of the standard Gauss algorithm @@ -291,8 +291,8 @@ END SUBROUTINE setdmc !> @param[in,out] b !> REAL(db), array, takes and returns the matrix b. !> @param[in] n -!> INTEGER, takes the number of grid points. -!--------------------------------------------------------------------------- +!> INTEGER, takes the number of grid points. +!--------------------------------------------------------------------------- PURE SUBROUTINE gauss(a,b,n) INTEGER,INTENT(IN) :: n REAL(db),INTENT(IN) :: a(n,n) @@ -313,7 +313,7 @@ PURE SUBROUTINE gauss(a,b,n) ! using d as intermediate storage IF(k/=ks) THEN d(k:2*n)=c(k,k:2*n) - c(k,k:2*n)=c(ks,k:2*n) + c(k,k:2*n)=c(ks,k:2*n) c(ks,k:2*n)=d(k:2*n) END IF c(k,k+1:2*n)=c(k,k+1:2*n)/c(k,k) diff --git a/Code/inout.f90 b/Code/inout.f90 index e4c199b..6e3d09c 100644 --- a/Code/inout.f90 +++ b/Code/inout.f90 @@ -1,12 +1,12 @@ !------------------------------------------------------------------------------ ! MODULE: Modulename !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief -!!This module contains the procedures for binary input and output -!!of the larger fields. +!!This module contains the procedures for binary input and output +!!of the larger fields. !> -!>@details +!>@details !!There are two variants, both of which are written !!at regular intervals: the wave function file \c wffile and the !!files containing the densities and currents. Since the former is @@ -42,13 +42,13 @@ MODULE Inout IMPLICIT NONE PRIVATE :: write_one_density,write_vec_density CONTAINS -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: write_wavefunctions !> @brief !!This subroutine writes the wave functions to the disk and has two modes of !!operation depending on whether the code runs on distributed-memory !!systems in \c MPI mode or on a shared-memory or single-processor -!!machine. +!!machine. !> !> @details !!In both cases it first determines the number of filled @@ -58,9 +58,9 @@ MODULE Inout !!states. !! !! - Sequential operation: -!! this case is recognized recognized by mpi_nprocs==1. Open +!! this case is recognized recognized by mpi_nprocs==1. Open !! \c wffile, then write four records containing general information. -!! - Record 1: iter, nstmax, nneut, nprot, number, npsi, +!! - Record 1: iter, nstmax, nneut, nprot, number, npsi, !! charge_number, mass_number, cm. !! - Record 2: nx, ny, nz, dx, dy, dz, wxyz. !! - Record 3: x, y, z. @@ -88,7 +88,7 @@ MODULE Inout !! processor and \c wffile in the form \c nnn.wffile. For !! example, if \c wffile has the value 'Ca40', these files will be !! \c 000.Ca40, \c 001.Ca40, \c 002.Ca40, etc. up to the number of processors. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE write_wavefunctions USE Parallel, ONLY: mpi_myproc,mpi_nprocs,nstloc,node,localindex INTEGER :: nst,iq,number(2) @@ -121,19 +121,19 @@ SUBROUTINE write_wavefunctions ENDDO CLOSE(UNIT=scratch2,STATUS='KEEP') END SUBROUTINE write_wavefunctions -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: write_densities !> @brief !!This subroutine produces a file \c iter.tdd with density and -!!current data for the present time step or iteration with number \c iter. -!!In the file name \c iter is given in 6 decimal digits. +!!current data for the present time step or iteration with number \c iter. +!!In the file name \c iter is given in 6 decimal digits. !> !> @details !!The record structure is as follows: !! -!! - Record 1: this contains the variables iter, time, nx, +!! - Record 1: this contains the variables iter, time, nx, !! ny, and nz to define the dimensions of the fields. -!! - Record 2: contains the variables dx, dy, dz, wxyz, x, +!! - Record 2: contains the variables dx, dy, dz, wxyz, x, !! y, and z to allow proper labelling of axes in plots, etc. !! - Further records: for each field to be written, a record is !! produced with the following information: @@ -148,8 +148,8 @@ END SUBROUTINE write_wavefunctions !! potential). !! . !! After this identification record, the corresponding field itself is -!! written. -!! The dimension varies in the following way: +!! written. +!! The dimension varies in the following way: !! !! !!
Dimensions of arrays
scalarwrite_isospindimension @@ -172,7 +172,7 @@ END SUBROUTINE write_wavefunctions !! -O: spin-orbit density \c sodens. Name \c s-o-Dens. !! !!This system is set up to be easily modified for writing additional fields. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE write_densities CHARACTER(10) :: filename CHARACTER(1) :: c @@ -205,7 +205,7 @@ SUBROUTINE write_densities END DO CLOSE(UNIT=scratch) END SUBROUTINE write_densities -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: write_one_density !> @brief !!This subroutines does the actual output for \c write_densities in @@ -215,7 +215,7 @@ END SUBROUTINE write_densities !> INTEGER, takes the name of the density. !> @param[in] values !> REAL(db), takes the density. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE write_one_density(name,values) CHARACTER(*),INTENT(IN) :: name REAL(db),INTENT(IN) :: values(nx,ny,nz,2) @@ -230,17 +230,17 @@ SUBROUTINE write_one_density(name,values) WRITE(scratch) a END IF END SUBROUTINE write_one_density -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: write_vec_density !> @brief -!!This also does the actual output for subroutines \c write_densities +!!This also does the actual output for subroutines \c write_densities !!for the case of a vector field. Its functioning should be clear from the !!description above. !> @param[in] name !> INTEGER, takes the name of the density. !> @param[in] values !> REAL(db), takes the vector density. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE write_vec_density(name,values) CHARACTER(*),INTENT(IN) :: name REAL(db),INTENT(IN) :: values(nx,ny,nz,3,2) @@ -255,7 +255,7 @@ SUBROUTINE write_vec_density(name,values) WRITE(scratch) a END IF END SUBROUTINE write_vec_density -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: plot_density !> @brief !!Produces a simple printer plot of the density distribution in the @@ -266,8 +266,8 @@ END SUBROUTINE write_vec_density !> @details !!It is based on a very old routine found at ORNL and was translated !!into modern Fortran. It uses helper function \c bplina for -!!interpolation. -!--------------------------------------------------------------------------- +!!interpolation. +!--------------------------------------------------------------------------- SUBROUTINE plot_density REAL(db),PARAMETER :: density_scale=0.14D0 INTEGER,PARAMETER :: ixsc=10,izsc=6 @@ -320,7 +320,7 @@ SUBROUTINE plot_density END IF WRITE(*,'(A,12(F6.2,4X),F6.2)') ' x= ',(xco(i),i=1,ntkx) END SUBROUTINE plot_density -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: bplina !> @brief !!Does a bilinear interpolation of fun(nx,nz), on xar(n) and zar(m) @@ -338,8 +338,8 @@ END SUBROUTINE plot_density !> @param[in] xcu !> REAL(db), takes the x-value at which the interpolation is performed !> @param[in] zcu -!> REAL(db), takes the z-value at which the interpolation is performed -!--------------------------------------------------------------------------- +!> REAL(db), takes the z-value at which the interpolation is performed +!--------------------------------------------------------------------------- PURE FUNCTION bplina(n,m,xar,zar,fun,xcu,zcu) RESULT(ff) INTEGER,INTENT(IN) :: n,m REAL(db),INTENT(IN) :: xar(n),zar(m),fun(n,m),xcu,zcu @@ -362,7 +362,7 @@ PURE FUNCTION bplina(n,m,xar,zar,fun,xcu,zcu) RESULT(ff) ff=dxf*(fun(icu+1,jcu+1)*dzf+fun(icu+1,jcu)*(1.0D0-dzf)) & + (1.D0-dxf)*(fun(icu,jcu+1)*dzf+fun(icu,jcu)*(1.0D0-dzf)) END FUNCTION bplina -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: sp_properties !> @brief !!In this routine the kinetic energy, orbital and spin angular momenta @@ -383,12 +383,12 @@ END FUNCTION bplina !!are combined to the desired matrix elements; the only technical point !!to remark is that since the result must be real, efficiency can be !!achieved by formulating the complex products in an explicit way. Then -!!\c kin contains the kinetic energy (without the \f$ \hbar^2/2m \f$), +!!\c kin contains the kinetic energy (without the \f$ \hbar^2/2m \f$), !!\c cc the orbital and \c ss then spin matrix elements. !! !!Finally only the volume element, the factor of one half for the -!!spin ad the prefactor of the kinetic energy are added. -!--------------------------------------------------------------------------- +!!spin ad the prefactor of the kinetic energy are added. +!--------------------------------------------------------------------------- SUBROUTINE sp_properties USE Trivial, ONLY: cmulx,cmuly,cmulz INTEGER :: nst,ix,iy,iz,is,ixx,iyy,izz @@ -397,7 +397,7 @@ SUBROUTINE sp_properties ALLOCATE(pst(nx,ny,nz,2),psx(nx,ny,nz,2),psy(nx,ny,nz,2),psz(nx,ny,nz,2),psw(nx,ny,nz,2)) sp_orbital=0.D0 sp_spin=0.D0 - sp_kinetic=0.0D0 + sp_kinetic=0.0D0 sp_parity=0.0D0 xx=x-cmtot(1) @@ -407,30 +407,30 @@ SUBROUTINE sp_properties IF(node(nst)/=mpi_myproc) CYCLE pst=psi(:,:,:,:,localindex(nst)) IF(TFFT) THEN - CALL cdervx(pst,psx) - CALL cdervy(pst,psy) + CALL cdervx(pst,psx) + CALL cdervy(pst,psy) - CALL cdervz(pst,psz) - CALL laplace(pst,psw) + CALL cdervz(pst,psz) + CALL laplace(pst,psw) ELSE - CALL cmulx(der1x,pst,psx,0) - CALL cmuly(der1y,pst,psy,0) - CALL cmulz(der1z,pst,psz,0) - CALL cmulx(der2x,pst,psw,0) - CALL cmuly(der2y,pst,psw,1) - CALL cmulz(der2z,pst,psw,1) + CALL cmulx(der1x,pst,psx,0) + CALL cmuly(der1y,pst,psy,0) + CALL cmulz(der1z,pst,psz,0) + CALL cmulx(der2x,pst,psw,0) + CALL cmuly(der2y,pst,psw,1) + CALL cmulz(der2z,pst,psw,1) ENDIF cc=0.D0 ss=0.D0 kin=0.D0 xpar=0.D0 - DO iz=1,nz - izz=nz-iz+1 - DO iy=1,ny - iyy=ny-iy+1 - DO ix=1,nx - ixx=nx-ix+1 - DO is=1,2 + DO iz=1,nz + izz=nz-iz+1 + DO iy=1,ny + iyy=ny-iy+1 + DO ix=1,nx + ixx=nx-ix+1 + DO is=1,2 rp=REAL(pst(ix,iy,iz,is)) ip=AIMAG(pst(ix,iy,iz,is)) cc(1)=cc(1)+ & @@ -462,19 +462,19 @@ SUBROUTINE sp_properties END DO DEALLOCATE(pst,psx,psy,psz,psw) END SUBROUTINE sp_properties -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: Routinename !> @brief !!This is given a file name and a character !!string for a header line to start the file contents. It is used for !!the *.res files. If the file already exists, nothing is done, !!since this probably a restart job and output should just be added at -!!the end of the file. +!!the end of the file. !> @param[in] filename !> CHARACTER, array, takes the filename. !> @param[in] header !> CHARACTER, array, takes the intended header for the file. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE start_protocol(filename,header) ! if the protocol file exists, do nothing, since later writes will be ! appended. Otherwise write the title and header lines into the new file diff --git a/Code/levels.f90 b/Code/levels.f90 index 0884d8c..7b17d42 100644 --- a/Code/levels.f90 +++ b/Code/levels.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! MODULE: Levels !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module is concerned with the wave function data: definition of !!the pertinent arrays, allocating their storage and simple operations on them. @@ -17,8 +17,8 @@ MODULE Levels !!wave functions are present on each node. Note that for all other !!wave-function related arrays, such as single-particle energies, the !!full set is stored on each node. - INTEGER :: nstloc ! npmin(1) through npsi(1) and the proton states !>run from npmin(2) through npsi(2). Protons follow - !>neutrons, so npmin(1)=1 and npmin(2)=npsi(1)+1. + !>neutrons, so npmin(1)=1 and npmin(2)=npsi(1)+1. !>Note that for each particle type the number of states can be !>larger than the particle number, as states may be fractionally !>occupied or even empty. If initialization is not from fragments, @@ -38,8 +38,8 @@ MODULE Levels !>meaning as the final index for proton states, which coincides with !>the total number of states, npsi(2)=nstmax. !>@{ - INTEGER :: npmin(2) - INTEGER :: npsi(2) + INTEGER :: npmin(2) + INTEGER :: npsi(2) !>@} REAL(db) :: charge_number,& ! (nstmax,nstmax) + COMPLEX(db), ALLOCATABLE :: hmatr(:,:) ! (nstmax,nstmax) !!and is used for the single-particle Hamiltonian in the diagonalization step. REAL(db), ALLOCATABLE, DIMENSION(:) :: sp_energy ! @brief !!This subroutine allocates all the arrays associated with -!!single-particle wave functions. +!!single-particle wave functions. !> !> @details !!Note that while most have dimension -!!\c nstmax, \c psi itself is dimensioned for the number \c nstloc of wave +!!\c nstmax, \c psi itself is dimensioned for the number \c nstloc of wave !!functions on one specific processor. It also records the isospin value. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE alloc_levels ALLOCATE(psi(nx,ny,nz,2,nstloc), & - sp_energy(nstmax),sp_efluct1(nstmax),sp_kinetic(nstmax),& + sp_energy(nstmax),sp_efluct1(nstmax),sp_kinetic(nstmax),& sp_norm(nstmax),sp_efluct2(nstmax),sp_parity(nstmax), & sp_orbital(3,nstmax),sp_spin(3,nstmax),wocc(nstmax), & isospin(nstmax)) isospin(1:npsi(1))=1 isospin(npmin(2):npsi(2))=2 END SUBROUTINE alloc_levels -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: cdervx !> @brief !!This routine calculates -!!derivatives of wave functions using the FFT method, in the \f$ x \f$-direction. +!!derivatives of wave functions using the FFT method, in the \f$ x \f$-direction. !> !> @details !!The last argument can be omitted @@ -111,11 +111,11 @@ END SUBROUTINE alloc_levels !> COMPLEX(db), array, returns the first derivative. !> @param[out] d2psout !> COMPLEX(db), array, OPTIONAL, returns the first derivative. -!--------------------------------------------------------------------------- - SUBROUTINE cdervx(psin,d1psout,d2psout) +!--------------------------------------------------------------------------- + SUBROUTINE cdervx(psin,d1psout,d2psout) COMPLEX(db), INTENT(IN) :: psin(:,:,:,:) COMPLEX(db), INTENT(OUT):: d1psout(:,:,:,:) - COMPLEX(db), INTENT(OUT), OPTIONAL :: d2psout(:,:,:,:) + COMPLEX(db), INTENT(OUT), OPTIONAL :: d2psout(:,:,:,:) REAL(db) :: kfac INTEGER :: ix kfac=(PI+PI)/(dx*nx) @@ -138,11 +138,11 @@ SUBROUTINE cdervx(psin,d1psout,d2psout) d1psout(nx/2+1,:,:,:)=(0.D0,0.D0) CALL dfftw_execute_dft(xbackward,d1psout,d1psout) END SUBROUTINE cdervx -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: cdervy !> @brief !!This routine calculates -!!derivatives of wave functions using the FFT method, in the \f$ y \f$-direction. +!!derivatives of wave functions using the FFT method, in the \f$ y \f$-direction. !> !> @details !!The last argument can be omitted @@ -156,11 +156,11 @@ END SUBROUTINE cdervx !> COMPLEX(db), array, returns the first derivative. !> @param[out] d2psout !> COMPLEX(db), array, OPTIONAL, returns the first derivative. -!--------------------------------------------------------------------------- - SUBROUTINE cdervy(psin,d1psout,d2psout) +!--------------------------------------------------------------------------- + SUBROUTINE cdervy(psin,d1psout,d2psout) COMPLEX(db), INTENT(IN) :: psin(:,:,:,:) COMPLEX(db), INTENT(OUT):: d1psout(:,:,:,:) - COMPLEX(db), INTENT(OUT), OPTIONAL :: d2psout(:,:,:,:) + COMPLEX(db), INTENT(OUT), OPTIONAL :: d2psout(:,:,:,:) REAL(db) :: kfac INTEGER :: iy,is,k kfac=(PI+PI)/(dy*ny) @@ -195,11 +195,11 @@ SUBROUTINE cdervy(psin,d1psout,d2psout) END DO END DO END SUBROUTINE cdervy -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: cdervz !> @brief !!This routine calculates -!!derivatives of wave functions using the FFT method, in the \f$ z \f$-direction. +!!derivatives of wave functions using the FFT method, in the \f$ z \f$-direction. !> !> @details !!The last argument can be omitted @@ -213,11 +213,11 @@ END SUBROUTINE cdervy !> COMPLEX(db), array, returns the first derivative. !> @param[out] d2psout !> COMPLEX(db), array, OPTIONAL, returns the first derivative. -!--------------------------------------------------------------------------- - SUBROUTINE cdervz(psin,d1psout,d2psout) +!--------------------------------------------------------------------------- + SUBROUTINE cdervz(psin,d1psout,d2psout) COMPLEX(db), INTENT(IN) :: psin(:,:,:,:) COMPLEX(db), INTENT(OUT):: d1psout(:,:,:,:) - COMPLEX(db), INTENT(OUT), OPTIONAL :: d2psout(:,:,:,:) + COMPLEX(db), INTENT(OUT), OPTIONAL :: d2psout(:,:,:,:) REAL(db) :: kfac INTEGER :: iz,is kfac=(PI+PI)/(dz*nz) @@ -246,7 +246,7 @@ SUBROUTINE cdervz(psin,d1psout,d2psout) CALL dfftw_execute_dft(zbackward,d1psout(:,:,:,is),d1psout(:,:,:,is)) END DO END SUBROUTINE cdervz -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: laplace !> @brief !!depending on the presence of the third argument \c e0inv, it can @@ -271,8 +271,8 @@ END SUBROUTINE cdervz !> COMPLEX(db), array, returns the resulting wave function. !> @param[in] e0inv !> REAL(db), OPTIONAL, takes the damping factor. -!--------------------------------------------------------------------------- - SUBROUTINE laplace(psin,psout,e0inv) +!--------------------------------------------------------------------------- + SUBROUTINE laplace(psin,psout,e0inv) USE Forces, ONLY: h2ma USE Grids, ONLY: dx,dy,dz COMPLEX(db), INTENT(IN) :: psin(:,:,:,:) @@ -317,19 +317,19 @@ SUBROUTINE laplace(psin,psout,e0inv) CALL dfftw_execute_dft(pbackward,psout(:,:,:,is),psout(:,:,:,is)) END DO END SUBROUTINE laplace -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: schmid !> @brief !!This performs a Gram-Schmidt orthogonalization of the single particle !!levels in the straightforward way: loop over states and subtract the -!!components of all lower-indexed states from a wave function. -!!the static calculations. +!!components of all lower-indexed states from a wave function. +!!the static calculations. !> !> @details !!It can be !!parallelized under \c OpenMP, albeit not very efficiently, but !!becomes too slow under \c MPI, so that \c MPI is not enabled for -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE schmid USE Trivial, ONLY: rpsnorm, overlap COMPLEX(db) :: oij,omptemp(nx,ny,nz,2) @@ -339,7 +339,7 @@ SUBROUTINE schmid IF(npsi(iq)-npmin(iq)+1 /= 1) THEN DO nst=npmin(iq),npsi(iq) omptemp=0 - !$OMP PARALLEL DO PRIVATE(j,oij) REDUCTION(+:omptemp) + !$OMP PARALLEL DO PRIVATE(j,oij) REDUCTION(+:omptemp) DO j=npmin(iq),nst-1 oij=overlap(psi(:,:,:,:,j),psi(:,:,:,:,nst)) omptemp(:,:,:,:)=omptemp(:,:,:,:)+oij*psi(:,:,:,:,j) diff --git a/Code/main3d.f90 b/Code/main3d.f90 index 305be6c..e7be85d 100644 --- a/Code/main3d.f90 +++ b/Code/main3d.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! Programm: tdhf3d !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This is the main program that organizes the reading of the input and !!funnels the calculation into the correct subroutines. @@ -9,7 +9,7 @@ !>@details !!It consists of a number of simple steps. They are: !! -# initialize \c MPI in case of an \c MPI parallel job. -!! Start reading from standard input beginning with +!! Start reading from standard input beginning with !! namelist files for any changes in the file names. !! -# read the definition of the force to be used and set !! it up. @@ -53,7 +53,7 @@ !! work. !! -# finally the \c MPI system is terminated. !------------------------------------------------------------------------------ -PROGRAM tdhf3d +PROGRAM tdhf3d USE Params USE Fourier USE Forces, ONLY: read_force @@ -173,7 +173,7 @@ PROGRAM tdhf3d CALL schmid WRITE(*,*) 'Reorthogonalization complete' END IF - ELSEIF(nof==0) THEN + ELSEIF(nof==0) THEN CALL harmosc ELSE CALL init_user diff --git a/Code/meanfield.f90 b/Code/meanfield.f90 index da0635e..85614eb 100644 --- a/Code/meanfield.f90 +++ b/Code/meanfield.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! MODULE: Meanfield !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module calculates all the ingredients needed for the !!energy functional and for applying the single-particle Hamiltonian to @@ -25,32 +25,32 @@ Module Meanfield USE Params, ONLY: db,tcoul USE Densities - USE Forces + USE Forces USE Grids, ONLY: nx,ny,nz,der1x,der2x,der1y,der2y,der1z,der2z USE Coulomb, ONLY: poisson,wcoul IMPLICIT NONE - REAL(db),ALLOCATABLE,DIMENSION(:,:,:,:) :: upot ! @brief !!This subroutine has the simple task of allocating all the fields that !!are local to the module \c Meanfield. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE alloc_fields ALLOCATE(upot(nx,ny,nz,2),bmass(nx,ny,nz,2),divaq(nx,ny,nz,2), & aq(nx,ny,nz,3,2),spot(nx,ny,nz,3,2),wlspot(nx,ny,nz,3,2), & @@ -62,7 +62,7 @@ SUBROUTINE alloc_fields divaq=0.D0 dbmass=0.D0 END SUBROUTINE alloc_fields -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: skyrme !> @brief !!In this subroutine the various fields are calculated from the @@ -83,7 +83,7 @@ END SUBROUTINE alloc_fields !!\f$ \longrightarrow \f$ (b1-b1p)*rho(:,:,:,iq)+b1*rho(:,:,:,ic) !!This decomposition is used in all applicable cases. !! -!!For intermediate results the fields \c workden (scalar) and \c workvec +!!For intermediate results the fields \c workden (scalar) and \c workvec !!(vector) are used. !! !!Now the subroutine proceeds in the following steps: @@ -100,14 +100,14 @@ END SUBROUTINE alloc_fields !! Slater exchange correction (only if the \c ex parameter in the !! force is nonzero) are added to \c upot for protons, iq=2 . !! -# the Laplacian is applied to the densities and -!! the result stored in \c workden. Then the remaining terms are constructed. +!! the result stored in \c workden. Then the remaining terms are constructed. !! Note that the \c iq -loop is combined with the following steps. !! -# the effective mass is calculated. !! -# the gradient of the density is calculated and the !! spin-orbit vector \f$ \vec W_q \f$ is constructed in \c wlspot. !! -# the curl of the spin density vector is calculated !! and stored in \c workvec. -!! -# the vector \f$ \vec A_q \f$ is calculated from the current density +!! -# the vector \f$ \vec A_q \f$ is calculated from the current density !! and the curl of the spin density. !! -# the curl of the current density is calculated and stored in \c spot. !! -# now the two isospin contributions in \c spot @@ -119,18 +119,18 @@ END SUBROUTINE alloc_fields !! \f$ B_q \f$ is calculated and stored in the vector variable \c dbmass. !! . !!This concludes the calculation of all scalar and vector fields needed -!!for the application of the Skyrme force. -!--------------------------------------------------------------------------- - SUBROUTINE skyrme +!!for the application of the Skyrme force. +!--------------------------------------------------------------------------- + SUBROUTINE skyrme USE Trivial, ONLY: rmulx,rmuly,rmulz - REAL(db),PARAMETER :: epsilon=1.0d-25 + REAL(db),PARAMETER :: epsilon=1.0d-25 REAL(db) :: rotspp,rotspn REAL(db),ALLOCATABLE :: workden(:,:,:,:),workvec(:,:,:,:,:) INTEGER :: ix,iy,iz,ic,iq,icomp ALLOCATE(workden(nx,ny,nz,2),workvec(nx,ny,nz,3,2)) ! Step 1: 3-body contribution to upot. - DO iq=1,2 - ic=3-iq + DO iq=1,2 + ic=3-iq upot(:,:,:,iq)=(rho(:,:,:,1)+rho(:,:,:,2))**f%power * & ((b3*(f%power+2.D0)/3.D0-2.D0*b3p/3.D0)*rho(:,:,:,iq) & +b3*(f%power+2.D0)/3.D0*rho(:,:,:,ic) & @@ -143,8 +143,8 @@ SUBROUTINE skyrme CALL rmuly(der1y,sodens(:,:,:,2,iq),workden(:,:,:,iq),1) CALL rmulz(der1z,sodens(:,:,:,3,iq),workden(:,:,:,iq),1) ENDDO - DO iq=1,2 - ic=3-iq + DO iq=1,2 + ic=3-iq upot(:,:,:,iq)=upot(:,:,:,iq) & -(b4+b4p)*workden(:,:,:,iq)-b4*workden(:,:,:,ic) ENDDO @@ -157,12 +157,12 @@ SUBROUTINE skyrme ENDIF ! Step 4: remaining terms of upot DO iq=1,2 - CALL rmulx(der2x,rho(:,:,:,iq),workden(:,:,:,iq),0) - CALL rmuly(der2y,rho(:,:,:,iq),workden(:,:,:,iq),1) + CALL rmulx(der2x,rho(:,:,:,iq),workden(:,:,:,iq),0) + CALL rmuly(der2y,rho(:,:,:,iq),workden(:,:,:,iq),1) CALL rmulz(der2z,rho(:,:,:,iq),workden(:,:,:,iq),1) ENDDO - DO iq=1,2 - ic=3-iq + DO iq=1,2 + ic=3-iq upot(:,:,:,iq)=upot(:,:,:,iq)+(b0-b0p)*rho(:,:,:,iq)+b0*rho(:,:,:,ic) & ! t1,t2, and tau-dependent part ! +(b1-b1p)*tau(:,:,:,iq)+b1*tau(:,:,:,ic) & @@ -181,7 +181,7 @@ SUBROUTINE skyrme (b4+b4p)*workvec(:,:,:,:,iq)+b4*workvec(:,:,:,:,ic) END DO ! Step 7: calculate curl of spin density vector, store in workvec - DO iq=1,2 + DO iq=1,2 CALL rmuly(der1y,sdens(:,:,:,3,iq),workvec(:,:,:,1,iq),0) CALL rmulz(der1z,sdens(:,:,:,2,iq),workvec(:,:,:,1,iq),-1) CALL rmulz(der1z,sdens(:,:,:,1,iq),workvec(:,:,:,2,iq),0) @@ -197,7 +197,7 @@ SUBROUTINE skyrme -(b4+b4p)*workvec(:,:,:,:,iq)-b4*workvec(:,:,:,:,ic) ENDDO ! Step 9: calculate the curl of the current density, stopr in spot - DO iq=1,2 + DO iq=1,2 CALL rmuly(der1y,current(:,:,:,3,iq),spot(:,:,:,1,iq),0) CALL rmulz(der1z,current(:,:,:,2,iq),spot(:,:,:,1,iq),-1) CALL rmulz(der1z,current(:,:,:,1,iq),spot(:,:,:,2,iq),0) @@ -206,40 +206,40 @@ SUBROUTINE skyrme CALL rmuly(der1y,current(:,:,:,1,iq),spot(:,:,:,3,iq),-1) ENDDO ! Step 10: combine isospin contributions - DO icomp=1,3 + DO icomp=1,3 DO iz=1,nz DO iy=1,ny - DO ix=1,nx - rotspp=spot(ix,iy,iz,icomp,1) - rotspn=spot(ix,iy,iz,icomp,2) - spot(ix,iy,iz,icomp,1)=-(b4+b4p)*rotspp-b4*rotspn - spot(ix,iy,iz,icomp,2)=-(b4+b4p)*rotspn-b4*rotspp + DO ix=1,nx + rotspp=spot(ix,iy,iz,icomp,1) + rotspn=spot(ix,iy,iz,icomp,2) + spot(ix,iy,iz,icomp,1)=-(b4+b4p)*rotspp-b4*rotspn + spot(ix,iy,iz,icomp,2)=-(b4+b4p)*rotspn-b4*rotspp ENDDO ENDDO ENDDO ENDDO - ! Step 11: calculate divergence of aq in divaq - DO iq=1,2 + ! Step 11: calculate divergence of aq in divaq + DO iq=1,2 CALL rmulx(der1x,aq(:,:,:,1,iq),divaq(:,:,:,iq),0) CALL rmuly(der1y,aq(:,:,:,2,iq),divaq(:,:,:,iq),1) CALL rmulz(der1z,aq(:,:,:,3,iq),divaq(:,:,:,iq),1) ENDDO ! Step 12: calculate the gradient of the effective mass in dbmass - DO iq=1,2 + DO iq=1,2 CALL rmulx(der1x,bmass(:,:,:,iq),dbmass(:,:,:,1,iq),0) CALL rmuly(der1y,bmass(:,:,:,iq),dbmass(:,:,:,2,iq),0) CALL rmulz(der1z,bmass(:,:,:,iq),dbmass(:,:,:,3,iq),0) ENDDO DEALLOCATE(workden,workvec) END SUBROUTINE skyrme -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: hpsi !> @brief !!This subroutine applies the single-particle Hamiltonian to a !!single-particle wave function \c pinn to produce an output wave !!function \c pout. The argument \c iq indicates the isospin for !!the wave function and \c eshift is an energy shift which is zero in -!!the dynamic calculation but crucial to the static algorithm (see +!!the dynamic calculation but crucial to the static algorithm (see !!\c grstep in module \c Static). !> !> @details @@ -274,18 +274,18 @@ END SUBROUTINE skyrme !! +\vec S_q\cdot\vec\sigma !! -\frac{\I}{2} \left[(\nabla\cdot\vec A_q)+2\vec A_q\cdot\nabla\right]. !!\f] -!! +!! !! -# the non-derivative parts not involving spin. These !! arise from \f$ U_q \f$ and \f$ -\tfrac{\I}{2}\,\nabla\cdot\vec A_q \f$, which !! are combined into a complex expression. The energy shift \c eshift is also included. !! -# the spin current coupling is constructed by simply !! using the explicit definition of the Pauli matrices and multiplying !! the resulting matrix onto the spinor wave function. -!! -# the first and second derivative in the \f$ x \f$-direction are evaluated -!! and stored in the arrays \c pswk and \c pswk2. The last term in the Hamiltonian +!! -# the first and second derivative in the \f$ x \f$-direction are evaluated +!! and stored in the arrays \c pswk and \c pswk2. The last term in the Hamiltonian !! gives rise to the two contributions -!! \f[ -\frac{\partial B_q}{\partial x}\frac{\partial}{\partial x}-B_q -!! \frac{\partial^2}{\partial x^2}, \f] +!! \f[ -\frac{\partial B_q}{\partial x}\frac{\partial}{\partial x}-B_q +!! \frac{\partial^2}{\partial x^2}, \f] !! of which the second is evaluated !! straightforwardly, while the first one is combined with the spin-orbit !! contribution. The part of \f$ \I\vec W_q\cdot(\vec\sigma\times\nabla) \f$ @@ -314,12 +314,12 @@ END SUBROUTINE skyrme !> COMPLEX(db), array, takes the wave function the Hamiltonian is supposed to be applied to. !> @param[out] pout !> COMPLEX(db), array, returns the Hamiltonian applied to the wave function. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE hpsi(iq,eshift,pinn,pout) USE Trivial, ONLY: cmulx, cmuly, cmulz USE Levels, ONLY: cdervx,cdervy,cdervz INTEGER :: iq - REAL(db) :: eshift + REAL(db) :: eshift COMPLEX(db),DIMENSION(:,:,:,:) :: pinn,pout INTENT(IN) :: iq,eshift INTENT(INOUT) :: pinn @@ -329,7 +329,7 @@ SUBROUTINE hpsi(iq,eshift,pinn,pout) COMPLEX(db),ALLOCATABLE,DIMENSION(:,:,:,:) :: pswk,pswk2 ALLOCATE(pswk(nx,ny,nz,2),pswk2(nx,ny,nz,2)) ! Step 1: non-derivative parts not involving spin - DO is=1,2 + DO is=1,2 pout(:,:,:,is)=CMPLX(upot(:,:,:,iq)-eshift, & -0D0,db)*pinn(:,:,:,is) ENDDO @@ -342,14 +342,14 @@ SUBROUTINE hpsi(iq,eshift,pinn,pout) *pinn(:,:,:,1) - spot(:,:,:,3,iq)*pinn(:,:,:,2) ! Step 3: derivative terms in x IF(TFFT) THEN - CALL cdervx(pinn,pswk,d2psout=pswk2) + CALL cdervx(pinn,pswk,d2psout=pswk2) ELSE - CALL cmulx(der1x,pinn,pswk,0) - CALL cmulx(der2x,pinn,pswk2,0) + CALL cmulx(der1x,pinn,pswk,0) + CALL cmulx(der2x,pinn,pswk2,0) ENDIF - DO is=1,2 - ic=3-is - sigis=(3-2*is)*0.5D0 + DO is=1,2 + ic=3-is + sigis=(3-2*is)*0.5D0 pout(:,:,:,is)=pout(:,:,:,is) & -CMPLX(dbmass(:,:,:,1,iq),0.5D0*aq(:,:,:,1,iq) & -sigis*wlspot(:,:,:,2,iq),db)*pswk(:,:,:,is) & @@ -363,21 +363,21 @@ SUBROUTINE hpsi(iq,eshift,pinn,pout) (aq(:,:,:,1,iq)+wlspot(:,:,:,2,iq))*pinn(:,:,:,2)& +0.5D0*wlspot(:,:,:,3,iq)*pinn(:,:,:,1) IF(TFFT) THEN - CALL cdervx(pswk2,pswk) + CALL cdervx(pswk2,pswk) ELSE - CALL cmulx(der1x,pswk2,pswk,0) + CALL cmulx(der1x,pswk2,pswk,0) ENDIF pout(:,:,:,:)=pout(:,:,:,:) + pswk(:,:,:,:) ! Step 4: derivative terms in y IF(TFFT) THEN - CALL cdervy(pinn,pswk,d2psout=pswk2) + CALL cdervy(pinn,pswk,d2psout=pswk2) ELSE - CALL cmuly(der1y,pinn,pswk,0) - CALL cmuly(der2y,pinn,pswk2,0) + CALL cmuly(der1y,pinn,pswk,0) + CALL cmuly(der2y,pinn,pswk2,0) ENDIF - DO is=1,2 - ic=3-is - sigis=(3-2*is)*0.5D0 + DO is=1,2 + ic=3-is + sigis=(3-2*is)*0.5D0 pout(:,:,:,is)=pout(:,:,:,is) & -CMPLX(dbmass(:,:,:,2,iq),0.5D0*aq(:,:,:,2,iq) & +sigis*wlspot(:,:,:,1,iq),db)*pswk(:,:,:,is) & @@ -391,21 +391,21 @@ SUBROUTINE hpsi(iq,eshift,pinn,pout) (aq(:,:,:,2,iq)-wlspot(:,:,:,1,iq))*pinn(:,:,:,2)& +CMPLX(0D0,0.5D0*wlspot(:,:,:,3,iq),db)*pinn(:,:,:,1) IF(TFFT) THEN - CALL cdervy(pswk2,pswk) + CALL cdervy(pswk2,pswk) ELSE - CALL cmuly(der1y,pswk2,pswk,0) + CALL cmuly(der1y,pswk2,pswk,0) ENDIF pout(:,:,:,:)=pout(:,:,:,:) + pswk(:,:,:,:) ! Step 5: derivative terms in z IF(TFFT) THEN - CALL cdervz(pinn,pswk,d2psout=pswk2) + CALL cdervz(pinn,pswk,d2psout=pswk2) ELSE - CALL cmulz(der1z,pinn,pswk,0) - CALL cmulz(der2z,pinn,pswk2,0) + CALL cmulz(der1z,pinn,pswk,0) + CALL cmulz(der2z,pinn,pswk2,0) ENDIF - DO is=1,2 - ic=3-is - sigis=(3-2*is)*0.5D0 + DO is=1,2 + ic=3-is + sigis=(3-2*is)*0.5D0 pout(:,:,:,is)=pout(:,:,:,is) & -CMPLX(dbmass(:,:,:,3,iq),0.5D0*aq(:,:,:,3,iq),db)*pswk(:,:,:,is) & +CMPLX(sigis*wlspot(:,:,:,1,iq),-0.5D0*wlspot(:,:,:,2,iq),db)* & @@ -416,9 +416,9 @@ SUBROUTINE hpsi(iq,eshift,pinn,pout) pswk2(:,:,:,2) = CMPLX(0D0,-0.5D0,db)*aq(:,:,:,3,iq)*pinn(:,:,:,2)& +CMPLX(-0.5D0*wlspot(:,:,:,1,iq),-0.5D0*wlspot(:,:,:,2,iq),db)*pinn(:,:,:,1) IF(TFFT) THEN - CALL cdervz(pswk2,pswk) + CALL cdervz(pswk2,pswk) ELSE - CALL cmulz(der1z,pswk2,pswk,0) + CALL cmulz(der1z,pswk2,pswk,0) ENDIF pout(:,:,:,:)=pout(:,:,:,:) + pswk(:,:,:,:) ! diff --git a/Code/moment.f90 b/Code/moment.f90 index cb2c8ea..e148686 100644 --- a/Code/moment.f90 +++ b/Code/moment.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! MODULE: Modulename !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!In this module various moments of the density distribution are !!calculated. Most of them come in two versions: an isospin-dependent @@ -11,7 +11,7 @@ !!stored: that of the neutrons cm(1:3,1) , of the protons !!cm(1:3,2), and of the total mass distribution cmtot(1:3) . !> -!>@details +!>@details !!Since the geometrical arrangement in space can be arbitrary with !!respect to the Cartesian coordinate system - this is certainly true !!for non-central collision situations - some quantities associated @@ -20,7 +20,7 @@ !!tensor and diagonalizes it to obtain the principal !!axes of the nucleus. In this frame then we compute the spherical !!quadrupole moments \f$ Q_{2m} \f$ with their -!!dimensionless counterparts \f$ a_o \f$, \f$ a_2 \f$ and the deformation +!!dimensionless counterparts \f$ a_o \f$, \f$ a_2 \f$ and the deformation !!parameters \f$ \beta \f$ and \f$ \gamma \f$. !------------------------------------------------------------------------------ @@ -31,11 +31,11 @@ MODULE Moment IMPLICIT NONE PRIVATE - REAL(db) :: pnr(2) !pnr(1)=\f$ N \f$ and + REAL(db) :: pnr(2) !pnr(1)=\f$ N \f$ and !!protons pnr(2)= \f$ Z \f$. These are obtained by a simple integration - !!of the densities \c rho. - REAL(db) :: pnrtot ! pnrtot=\f$ A \f$. - !!These are obtained by a simple integration of the densities \c rho. + !!of the densities \c rho. + REAL(db) :: pnrtot ! pnrtot=\f$ A \f$. + !!These are obtained by a simple integration of the densities \c rho. REAL(db) :: cm(3,2) ! @brief !!This is the principal subroutine for calculating the geometric @@ -112,15 +112,15 @@ MODULE Moment !! components \c qmat of the Cartesian quadrupole tensor !! \f$ Q_{kl}=\int \D^3r \left(3x_kx_l-r^2\delta_{kl}\right) \rho(\vec r) \f$. !! - After this, the subroutine \c q2diag is used to determine the -!! spherical components \f$ Q_{20} \f$ and \f$ Q_{22} \f$ in the +!! spherical components \f$ Q_{20} \f$ and \f$ Q_{22} \f$ in the !! principal-axes frame also generating some printout in the process. !! These are then multiplied with a scale factor to yield the -!! dimensionless deformation parameters \f$ a_0 \f$ and \f$ a_2 \f$ and -!! finally by conversion to polar coordinates the Bohr-Mottelson parameters +!! dimensionless deformation parameters \f$ a_0 \f$ and \f$ a_2 \f$ and +!! finally by conversion to polar coordinates the Bohr-Mottelson parameters !! \f$ \beta \f$ and \f$ \gamma \f$. !! - The Cartesian and polar deformation parameters are then printed. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE moments(L_val,M_val) USE Densities, ONLY: rho,current INTEGER :: ix,iy,iz,iq,L_val,M_val,iii @@ -131,12 +131,12 @@ SUBROUTINE moments(L_val,M_val) pnr=0.D0 cm=0.D0 pcm=0.D0 - DO iq=1,2 - DO iz=1,nz - xx(3)=z(iz) - DO iy=1,ny - xx(2)=y(iy) - DO ix=1,nx + DO iq=1,2 + DO iz=1,nz + xx(3)=z(iz) + DO iy=1,ny + xx(2)=y(iy) + DO ix=1,nx xx(1)=x(ix) pnr(iq)=pnr(iq)+wxyz*rho(ix,iy,iz,iq) cm(:,iq)=cm(:,iq)+wxyz*xx*rho(ix,iy,iz,iq) @@ -145,8 +145,8 @@ SUBROUTINE moments(L_val,M_val) ENDDO ENDDO ENDDO - pnrtot=pnr(1)+pnr(2) - cmtot=(cm(:,1)+cm(:,2))/pnrtot + pnrtot=pnr(1)+pnr(2) + cmtot=(cm(:,1)+cm(:,2))/pnrtot DO iq=1,2 cm(:,iq)=cm(:,iq)/pnr(iq) ENDDO @@ -165,17 +165,17 @@ SUBROUTINE moments(L_val,M_val) HexaDeca = 0.0D0 DiaTriaConta = 0.0D0 - DO iq=1,2 - DO iz=1,nz - xx(3)=z(iz)-cm(3,iq) - x2(3)=xx(3)**2 - DO iy=1,ny - xx(2)=y(iy)-cm(2,iq) + DO iq=1,2 + DO iz=1,nz + xx(3)=z(iz)-cm(3,iq) + x2(3)=xx(3)**2 + DO iy=1,ny + xx(2)=y(iy)-cm(2,iq) x2(2)=xx(2)**2 - DO ix=1,nx - xx(1)=x(ix)-cm(1,iq) - x2(1)=xx(1)**2 - vol=wxyz*rho(ix,iy,iz,iq) + DO ix=1,nx + xx(1)=x(ix)-cm(1,iq) + x2(1)=xx(1)**2 + vol=wxyz*rho(ix,iy,iz,iq) rms(iq)=vol*SUM(x2)+rms(iq) r3(iq)=vol*(SUM(x2)**(1.5d0))+r3(iq) r4(iq)=vol*(SUM(x2)**(2))+r4(iq) @@ -186,22 +186,22 @@ SUBROUTINE moments(L_val,M_val) qmat(2,2,iq)=qmat(2,2,iq)+vol*(x2(2)+x2(2)-x2(1)-x2(3)) qmat(2,3,iq)=qmat(2,3,iq)+3.D0*vol*xx(2)*xx(3) qmat(3,3,iq)=qmat(3,3,iq)+vol*(x2(3)+x2(3)-x2(1)-x2(2)) - x2m(:,iq)=vol*x2(:)+x2m(:,iq) + x2m(:,iq)=vol*x2(:)+x2m(:,iq) Mono(iq)=vol*((0.5d0*SQRT(1.0d0/PI))*(SUM(x2)))+Mono(iq) - + eta = vol*(SUM(x2))*5.0d0/3.0d0 Di_is(iq)=vol*(Y_lm(1,M_val,x(ix),y(iy),z(iz))*(SQRT(x(ix)**2+y(iy)**2+z(iz)**2)**3- & SQRT(x(ix)**2+y(iy)**2+z(iz)**2)*eta)*SQRT(2*1+1.0d0))+Di_is(iq) ! print*,'inside loop 1' - + tmp=(Y_lm(1,M_val,x(ix),y(iy),z(iz))*SQRT(3.0d0)) tmp = tmp*SQRT(x(ix)**2+y(iy)**2+z(iz)**2)*vol Di_iv(iq)=tmp+Di_iv(iq) ! print*,Di(iq),tmp,iq,'inside looop 2' - + Quad(iq)=vol*(Y_lm(2,M_val,xx(1),xx(2),xx(3))*SUM(x2)*SQRT(2*2+1.0d0))+Quad(iq) Oct(iq)=vol*(Y_lm(3,M_val,xx(1),xx(2),xx(3))*(SQRT(SUM(x2))**3) & @@ -218,21 +218,21 @@ SUBROUTINE moments(L_val,M_val) qmat(3,2,iq)=qmat(2,3,iq) ! print*,Di(iq),iq ENDDO - + r2tot=(rms(1)+rms(2))/pnrtot rmstot=SQRT((rms(1)+rms(2))/pnrtot) - rms=SQRT(rms/pnr) - x2mtot=(x2m(:,1)+x2m(:,2))/pnrtot + rms=SQRT(rms/pnr) + x2mtot=(x2m(:,1)+x2m(:,2))/pnrtot r3tot=(r3(1)+r3(2))/pnrtot r4tot=(r4(1)+r4(2))/pnrtot do iii=1,2 - IF(iii==1) THEN - facn=1.0D0 - facp=1.0D0 + IF(iii==1) THEN + facn=1.0D0 + facp=1.0D0 Di_tot(iii) = facn*Di_is(1)+facp*Di_is(2) - ELSE + ELSE if (L_val .eq. 1)then facn=-(pnr(2)/pnrtot) facp=(pnr(1))/pnrtot @@ -246,13 +246,13 @@ SUBROUTINE moments(L_val,M_val) Quad_tot(iii) = facn*Quad(1)+facp*Quad(2) Oct_tot(iii) = facn*Oct(1)+facp*Oct(2) HexDec_tot(iii) = facn*HexaDeca(1)+facp*HexaDeca(2) - DiaTriaConta_tot(iii) = facn*DiaTriaConta(1)+facp*DiaTriaConta(2) + DiaTriaConta_tot(iii) = facn*DiaTriaConta(1)+facp*DiaTriaConta(2) end do DO iq=1,2 x2m(:,iq)=x2m(:,iq)/pnr(iq) ENDDO - qmtot=qmat(:,:,1)+qmat(:,:,2) + qmtot=qmat(:,:,1)+qmat(:,:,2) IF(printnow.AND.wflag) WRITE(*,'(/A)') 'Cartesian quadrupole tensor,& & principal values, and axes:' CALL q2diag(qmat(:,:,1),q20(1),q22(1),'Neutrons ') @@ -272,54 +272,54 @@ SUBROUTINE moments(L_val,M_val) ' Beta20: ',beta20tot,' Beta22: ',beta22tot,' Beta: ',beta, & ' Gamma: ',gamma END SUBROUTINE moments -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: moment_shortprint !> @brief !!This subroutine simply prints some information into the specialized !!output files. \c monopolesfile receives the monopole moment, while \c quadrupolesfile !!receives the quadrupole moment and so on for other multipole files. The physical time !!starts each line to enable easy time-curve plotting. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE moment_shortprint() - OPEN(unit=scratch,file=monopolesfile,POSITION='APPEND') + OPEN(unit=scratch,file=monopolesfile,POSITION='APPEND') WRITE(scratch,'(F25.5,2x,F25.15,2x,F25.15)') time,Mono_tot(1),Mono_tot(2) CLOSE(unit=scratch) - - OPEN(unit=scratch,file=dipolesfile,POSITION='APPEND') + + OPEN(unit=scratch,file=dipolesfile,POSITION='APPEND') WRITE(scratch,'(F25.5,2x,F25.15,2x,F25.15)') time,Di_tot(1),Di_tot(2) CLOSE(unit=scratch) - OPEN(unit=scratch,file=quadrupolesfile,POSITION='APPEND') + OPEN(unit=scratch,file=quadrupolesfile,POSITION='APPEND') WRITE(scratch,'(F25.5,2x,F25.15,2x,F25.15)') time,Quad_tot(1),Quad_tot(2) CLOSE(unit=scratch) - - OPEN(unit=scratch,file=octupolesfile,POSITION='APPEND') + + OPEN(unit=scratch,file=octupolesfile,POSITION='APPEND') WRITE(scratch,'(F25.5,2x,F25.15,2x,F25.15)') time,Oct_tot(1),Oct_tot(2) CLOSE(unit=scratch) - OPEN(unit=scratch,file=hexadecapolesfile,POSITION='APPEND') + OPEN(unit=scratch,file=hexadecapolesfile,POSITION='APPEND') WRITE(scratch,'(F25.5,2x,F25.15,2x,F25.15)') time,HexDec_tot(1),HexDec_tot(2) CLOSE(unit=scratch) - OPEN(unit=scratch,file=diatriacontapolesfile,POSITION='APPEND') + OPEN(unit=scratch,file=diatriacontapolesfile,POSITION='APPEND') WRITE(scratch,'(F25.5,2x,F25.15,2x,F15.10)') time,DiaTriaConta_tot(1),DiaTriaConta_tot(2) CLOSE(unit=scratch) END SUBROUTINE moment_shortprint -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: moment_print !> @brief !!This subroutine prints a somewhat more detailed information. Particle !!number , r.m.s. radius, \f$ Q_{20} \f$, \c x2m, and center-of-mass are !!printed for the total distribution and also separately for neutrons !!and protons. This output goes to the regular output unit. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE moment_print INTEGER :: iq CHARACTER(11),PARAMETER :: Name(2)=(/ ' Neutron: ',' Proton: '/) Write(*,'(A)') ' Part.Num. rms-radius q20 & & & - & ' + & ' WRITE(*,'(a,2f12.4,1p,4e12.4,3e15.7)') ' Total: ',pnrtot,rmstot, & q20tot,x2mtot,cmtot DO iq=1,2 @@ -327,7 +327,7 @@ SUBROUTINE moment_print x2m(:,iq),cm(:,iq) ENDDO END SUBROUTINE moment_print -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: q2diag !> @brief !!his subroutine diagonalizes the Cartesian quadrupole tensor. To this @@ -341,17 +341,17 @@ END SUBROUTINE moment_print !!quadrupole moment: !!\f[ Q_{20}=\sqrt{\frac{5}{16\pi}}\,Q_{zz},\qquad !!Q_{22}=\sqrt{\frac{5}{96\pi}}\,\left(Q_{yy}-Q_{xx}\right). \f] They are -!!returned in \c q20x and \c q22x. +!!returned in \c q20x and \c q22x. !> !> @param[in,out] q_mat !> REAL(db), array, takes the quadrupole matrix. !> @param[out] q20x !> REAL(db), returns the value \f$ Q_{20} \f$. !> @param[out] q22x -!> REAL(db), returns the value \f$ Q_{22} \f$. +!> REAL(db), returns the value \f$ Q_{22} \f$. !> @param[in] title -!> CHARACTER, array, takes a title for printout. -!--------------------------------------------------------------------------- +!> CHARACTER, array, takes a title for printout. +!--------------------------------------------------------------------------- SUBROUTINE q2diag(q_mat,q20x,q22x,title) REAL(db),INTENT(INOUT) :: q_mat(3,3) REAL(db),INTENT(OUT) :: q20x,q22x @@ -359,7 +359,7 @@ SUBROUTINE q2diag(q_mat,q20x,q22x,title) REAL(db) :: q_eig(3),q_vec(3,3),fv1(20) INTEGER :: info,i, j,k - if(printnow.AND.wflag) write(*,'(3(f12.5,1x))') ((q_mat(j,k),k=1,3),j=1,3) + if(printnow.AND.wflag) write(*,'(3(f12.5,1x))') ((q_mat(j,k),k=1,3),j=1,3) CALL DSYEV('V','U',3,q_mat,3,q_eig,fv1,20,info) q_vec=q_mat IF(info/=0) STOP 'Quadrupole diagonalization failed' diff --git a/Code/pairs.f90 b/Code/pairs.f90 index 7b3b748..3f5f475 100644 --- a/Code/pairs.f90 +++ b/Code/pairs.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! MODULE: Pairs !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!The principal part of this module is the subroutine \c pair, which !!computes the pairing solution based on the BCS model. It is the only @@ -38,7 +38,7 @@ MODULE Pairs REAL(db),SAVE :: epair(2) ! @brief !!This is the only routine visible from outside the module. It solves @@ -61,32 +61,32 @@ MODULE Pairs !!two isospin values \c pairdn is called with the correct particle !!number as argument. This does the real work of solving the equations. !!Finally summary information is printed. -!--------------------------------------------------------------------------- - SUBROUTINE pair +!--------------------------------------------------------------------------- + SUBROUTINE pair REAL(db) :: particle_number ! prepare phase space weight for subroutines IF(.NOT.ALLOCATED(deltaf)) ALLOCATE(deltaf(nstmax)) ! calculate gaps CALL pairgap ! solve pairing problem for each isospin iq - DO iq=1,2 + DO iq=1,2 IF(iq==2) THEN particle_number=charge_number - ELSE + ELSE particle_number=mass_number-charge_number ENDIF CALL pairdn(particle_number) ENDDO ! print pairing information - IF(printnow.AND.wflag) THEN + IF(printnow.AND.wflag) THEN WRITE(*,'(/7x,a)') ' e_ferm e_pair aver_gap aver_force ' - DO iq=1,2 + DO iq=1,2 WRITE(*,'(a,i2,a,4(1pg12.4))') 'iq=',iq,': ',eferm(iq) , & epair(iq) ,avdelt(iq),avg(iq) ENDDO ENDIF END SUBROUTINE pair -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: pairgap !> @brief !!This subroutine calculates the pairing gaps \f$ \Delta_k \f$ stored in the @@ -98,7 +98,7 @@ END SUBROUTINE pair !!calculation. All gaps are set equal to \f$ 11.2\,{\rm MeV}/\sqrt{A} \f$ in !!this case. !! -!!In the general case, there is a loop over the two isospin values +!!In the general case, there is a loop over the two isospin values !!\c iq. The pairing density is obtained by evaluation of !!\f[ {\tt work}(\vec r)=\sum_k u_kv_k\left|\phi_k(\vec r)\right|^2 \f] !!where the simple conversion @@ -111,62 +111,62 @@ END SUBROUTINE pair !!the pairing density is simply multiplied by \c v0act, while for !!\c DDDI pairing ( ipair=6) it is !!\f[ V_P(\vec r)={\tt v0act}\cdot{\tt work}\cdot (1-x\rho(\vec r)/{\tt rho0pr}) \f] -!!involving the total density \f$ \rho \f$ and the parameter -!!\f[ rho0pr \f] from the pairing force definition and the parameter \f$x\f$ is used to create -!! mix pairing when its value is set to 0.5, it is read from the pairing definition in the +!!involving the total density \f$ \rho \f$ and the parameter +!!\f[ rho0pr \f] from the pairing force definition and the parameter \f$x\f$ is used to create +!! mix pairing when its value is set to 0.5, it is read from the pairing definition in the !! forces.data file. !! !!In the final step the gaps are computed as the expectation values of !!the pairing field, !!\f[ \Delta_k=\int\D^3r V_P(\vec r)\left|\phi_k(\vec r)\right|^2. \f] -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE pairgap ! computes the pair density INTEGER,PARAMETER :: itrsin=10 ! iqq is different from module variable iq INTEGER :: iqq,nst,is - REAL(db),PARAMETER :: smallp=0.000001 + REAL(db),PARAMETER :: smallp=0.000001 REAL(db) :: v0act - REAL(db) :: work(nx,ny,nz) + REAL(db) :: work(nx,ny,nz) ! constant gap in early stages of iteration - IF(iter<=itrsin) THEN + IF(iter<=itrsin) THEN deltaf=11.2/SQRT(mass_number) - RETURN + RETURN ENDIF ! now the detailed gaps: - DO iqq=1,2 + DO iqq=1,2 work=0.D0 ! accumulate new pair-density - DO nst=npmin(iqq),npsi(iqq) - DO is=1,2 + DO nst=npmin(iqq),npsi(iqq) + DO is=1,2 work=work+SQRT(MAX(wocc(nst)-wocc(nst)**2,smallp))*0.5* & (REAL(psi(:,:,:,is,nst))**2+AIMAG(psi(:,:,:,is,nst))**2) ENDDO ENDDO ! determine pairing strength - IF(iqq==2) THEN - v0act=p%v0prot - ELSE - v0act=p%v0neut + IF(iqq==2) THEN + v0act=p%v0prot + ELSE + v0act=p%v0neut ENDIF ! now multiply with strength to obtain local pair-potential IF(ipair==6) THEN work=v0act*work*(1D0-p%mixture*(rho(:,:,:,1)+rho(:,:,:,2))/p%rho0pr) ELSE - work=v0act*work + work=v0act*work END IF ! finally compute the actual gaps as s.p. expectation values with ! the pair potential - DO nst=npmin(iqq),npsi(iqq) - deltaf(nst)=0.0D0 - DO is=1,2 + DO nst=npmin(iqq),npsi(iqq) + deltaf(nst)=0.0D0 + DO is=1,2 deltaf(nst)=deltaf(nst)+wxyz*SUM(work* & (REAL(psi(:,:,:,is,nst))**2+AIMAG(psi(:,:,:,is,nst))**2)) ENDDO ENDDO ENDDO END SUBROUTINE pairgap -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: Routinename !> @brief !!The subroutine \c pairdn determines the pairing solution by using @@ -180,41 +180,41 @@ END SUBROUTINE pairgap !!single-particle energies. Then \c rbrent is called to calculate the !!correct solution, after which there is a loop for a straightforward !!evaluation of the module variables \c epair, \c avdelt, -!!and \c avg. +!!and \c avg. !> @param[in] particle_number -!> REAL(db), takes the particle number. -!--------------------------------------------------------------------------- +!> REAL(db), takes the particle number. +!--------------------------------------------------------------------------- SUBROUTINE pairdn(particle_number) REAL(db),INTENT(IN) :: particle_number REAL(db),PARAMETER :: xsmall=1.d-20 INTEGER :: it,na REAL(db) :: sumuv,sumduv,sumepa,edif,equasi,v2,vol ! start with non-pairing value of Fermi energy - it=npmin(iq)+NINT(particle_number)-1 - eferm(iq)=0.5D0*(sp_energy(it)+sp_energy(it+1)) + it=npmin(iq)+NINT(particle_number)-1 + eferm(iq)=0.5D0*(sp_energy(it)+sp_energy(it+1)) wocc(npmin(iq):it)=1.D0 wocc(it+1:npsi(iq))=0.D0 ! determine true fermi energy for given 'delta' and 'e' eferm(iq)=rbrent(particle_number) ! finally compute average gap and pairing energy - sumuv=0.0D0 - sumduv=0.0D0 + sumuv=0.0D0 + sumduv=0.0D0 sumepa=0.0D0 DO na=npmin(iq),npsi(iq) - edif=sp_energy(na)-eferm(iq) - equasi=SQRT(edif*edif+deltaf(na)**2) - v2=0.5D0-0.5D0*edif/equasi - vol=0.5D0*SQRT(MAX(v2-v2*v2,xsmall)) - sumuv=vol+sumuv - sumduv=vol*deltaf(na)+sumduv - sumepa=0.5D0*deltaf(na)**2/equasi+sumepa + edif=sp_energy(na)-eferm(iq) + equasi=SQRT(edif*edif+deltaf(na)**2) + v2=0.5D0-0.5D0*edif/equasi + vol=0.5D0*SQRT(MAX(v2-v2*v2,xsmall)) + sumuv=vol+sumuv + sumduv=vol*deltaf(na)+sumduv + sumepa=0.5D0*deltaf(na)**2/equasi+sumepa ENDDO - sumuv=MAX(sumuv,xsmall) - avdelt(iq)=sumduv/sumuv - epair(iq)=sumduv - avg(iq)=epair(iq)/sumuv**2 + sumuv=MAX(sumuv,xsmall) + avdelt(iq)=sumduv/sumuv + epair(iq)=sumduv + avg(iq)=epair(iq)/sumuv**2 END SUBROUTINE pairdn -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: rbrent !> @brief !!This subroutine is an adapted version of the Van @@ -228,8 +228,8 @@ END SUBROUTINE pairdn !!Fortran; it will be replaced at some time in the future. !> !> @param[in] particle_number -!> REAL(db), takes the particle number. -!--------------------------------------------------------------------------- +!> REAL(db), takes the particle number. +!--------------------------------------------------------------------------- FUNCTION rbrent(particle_number) RESULT(res) REAL(db),INTENT(in) :: particle_number REAL(db) :: res @@ -238,105 +238,105 @@ FUNCTION rbrent(particle_number) RESULT(res) REAL(db) :: ra,rb,rfa,rfb,rfc,rc,rd,re,rtol1,rxm,rp,rq,rr,rs, & rmin1,rmin2 INTEGER :: ii - ra=ra0 + ra=ra0 rb=rb0 CALL bcs_occupation(ra,rfa) rfa=particle_number-rfa CALL bcs_occupation(rb,rfb) rfb=particle_number-rfb - rfc=rfb - rc=rb - rd=rb-ra - re=rd - IF((rfa>0.D0) .AND.(rfb>0.D0)) THEN - IF(rfa>rfb) THEN - rb=1.0 - ELSE - rb=0.D0 + rfc=rfb + rc=rb + rd=rb-ra + re=rd + IF((rfa>0.D0) .AND.(rfb>0.D0)) THEN + IF(rfa>rfb) THEN + rb=1.0 + ELSE + rb=0.D0 ENDIF - ELSEIF((rfa<0.D0) .AND.(rfb<0.D0)) THEN - IF(rfa>rfb) THEN - rb=0.D0 - ELSE - rb=1.0 + ELSEIF((rfa<0.D0) .AND.(rfb<0.D0)) THEN + IF(rfa>rfb) THEN + rb=0.D0 + ELSE + rb=1.0 ENDIF - ELSE - iteration: DO ii=1,iitmax + ELSE + iteration: DO ii=1,iitmax IF(((rfb>0.D0) .AND.(rfc>0.D0)).OR.((rfb<0.D0) & .AND.(rfc<0.D0))) THEN - rc=ra - rfc=rfa - rd=rb-ra - re=rd + rc=ra + rfc=rfa + rd=rb-ra + re=rd ENDIF - IF(ABS(rfc)=rtol1).OR.(ABS(rfa)>ABS(rfb))) & THEN - rs=rfb/rfa - IF(ra==rc) THEN - rp=2.0*rxm*rs - rq=1.0-rs - ELSE - rq=rfa/rfc - rr=rfb/rfc + rs=rfb/rfa + IF(ra==rc) THEN + rp=2.0*rxm*rs + rq=1.0-rs + ELSE + rq=rfa/rfc + rr=rfb/rfc rp=rs *(2.0*rxm*rq *(rq-rr)-(rb-ra) *(rr- & 1.0)) - rq=(rq-1.0) *(rr-1.0) *(rs-1.0) + rq=(rq-1.0) *(rr-1.0) *(rs-1.0) ENDIF - IF(rp>0.D0) THEN - rq=-rq + IF(rp>0.D0) THEN + rq=-rq ENDIF - rp=ABS(rp) - rmin1=3.0*rxm*rq-ABS(rtol1*rq) - rmin2=ABS(rq*re) - IF(2.0*rprtol1) THEN - rb=rb+rd - ELSE - rb=rb+SIGN(rtol1,rxm) + ra=rb + rfa=rfb + IF(ABS(rd) >rtol1) THEN + rb=rb+rd + ELSE + rb=rb+SIGN(rtol1,rxm) ENDIF CALL bcs_occupation(rb,rfb) rfb=particle_number-rfb ENDDO iteration - STOP 'No solution found in pairing iterations' + STOP 'No solution found in pairing iterations' ENDIF END FUNCTION rbrent -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: Routinename !> @brief !!For a given Fermi energy \f$ \epsilon_F \f$ passed as argument \c efermi, !!this subroutine evaluates the particle number that would result with -!!such a Fermi energy and returns it as its second argument, \c bcs_partnum. +!!such a Fermi energy and returns it as its second argument, \c bcs_partnum. !> !> @details -!!The isospin is controlled by module variable \c iq. First the occupation +!!The isospin is controlled by module variable \c iq. First the occupation !!probabilities are calculated using the standard BCS expression !!\f[ v_k^2=\frac1{2}\left(1-\frac{\epsilon_k-\epsilon_F} !! {\sqrt{(\epsilon_k-\epsilon_F)^2+\Delta_k^2}}\right). \f] @@ -348,18 +348,18 @@ END FUNCTION rbrent !> REALD(db), takes the fermi energy. !> @param[out] bcs_partnum !> REAL(db), returns the oarticle number -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE bcs_occupation(efermi,bcs_partnum) - REAL(db),PARAMETER :: smal=1.0d-10 + REAL(db),PARAMETER :: smal=1.0d-10 REAL(db),INTENT(in) :: efermi REAL(db),INTENT(OUT) :: bcs_partnum INTEGER :: k REAL(db) :: edif - bcs_partnum=0.0 - DO k=npmin(iq),npsi(iq) - edif=sp_energy(k)-efermi - wocc(k)=0.5D0 *(1.0D0-edif/SQRT(edif**2+deltaf(k)**2)) - wocc(k)=MIN(MAX(wocc(k),smal),1.D0-smal) + bcs_partnum=0.0 + DO k=npmin(iq),npsi(iq) + edif=sp_energy(k)-efermi + wocc(k)=0.5D0 *(1.0D0-edif/SQRT(edif**2+deltaf(k)**2)) + wocc(k)=MIN(MAX(wocc(k),smal),1.D0-smal) bcs_partnum=bcs_partnum+wocc(k) ENDDO END SUBROUTINE bcs_occupation diff --git a/Code/parallel.f90 b/Code/parallel.f90 index e552934..96b39d8 100644 --- a/Code/parallel.f90 +++ b/Code/parallel.f90 @@ -1,14 +1,14 @@ !------------------------------------------------------------------------------ ! MODULE: Parallel !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module organizes the execution on distributed-memory machines !!using \c MPI. Its interface is made such that on sequential !!machines \c parallel.f90 can simply be replaced by a special !!version \c sequential.f90 which contains a definition of this !!module generating trivial data. -!> +!> !> @details !!\c MPI parallelization is !!based on distributing the wave functions onto the different nodes. @@ -27,8 +27,8 @@ !!communicate many small arrays. !! !!Since an efficient way of dealing with Gram-Schmidt orthogonalization -!!in this case was yet not found, at present the code can be run in \c MPI -!!parallel mode only for the dynamic case. +!!in this case was yet not found, at present the code can be run in \c MPI +!!parallel mode only for the dynamic case. !------------------------------------------------------------------------------ MODULE Parallel USE Params, ONLY : wflag @@ -47,30 +47,30 @@ MODULE Parallel !!\c node(i) and its index on that node is \c localindex(i). INTEGER, ALLOCATABLE :: globalindex(:) ! on the local node, i=1..nstloc, the - !!single-particle energy must be obtained using sp_energy(globalindex(i)). + !!single-particle energy must be obtained using sp_energy(globalindex(i)). INTEGER :: mpi_nprocs ! @brief This subroutine merely allocates the internal arrays of module \c Parallel. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE alloc_nodes ALLOCATE(node(nstmax),localindex(nstmax),globalindex(nstmax)) END SUBROUTINE alloc_nodes -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: init_all_mpi !> @brief !!\c init_all_mpi This subroutine initializes \c MPI and !!finds out the number of processors \c mpi_nprocs as well as the !!index of the current one \c mpi_myproc. The flag \c wflag is !!set to true only for the processor numbered 0. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE init_all_mpi CALL mpi_init(mpi_ierror) CALL mpi_comm_size(mpi_comm_world,mpi_nprocs,mpi_ierror) @@ -79,7 +79,7 @@ SUBROUTINE init_all_mpi IF(wflag) WRITE(*,'(a,i5)') ' number of nodes=',mpi_nprocs CALL mpi_barrier (mpi_comm_world, mpi_ierror) END SUBROUTINE init_all_mpi -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: associate_nodes !> @brief !!The first loop in this subroutine distributes the wave functions over @@ -93,11 +93,11 @@ END SUBROUTINE init_all_mpi !!complete sequence. The third loop sets up the reverse pointers !!\c localindex, which has to be done in a loop over all processors to !!set up the information for the proper global indices. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE associate_nodes INTEGER :: ncount,nst,ip,iloc ncount=0 - DO nst=1,nstmax + DO nst=1,nstmax ncount=MOD(ncount,mpi_nprocs) node(nst)=ncount ncount=ncount+1 @@ -124,14 +124,14 @@ SUBROUTINE associate_nodes ENDIF CALL mpi_barrier (mpi_comm_world, mpi_ierror) END SUBROUTINE associate_nodes -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: collect_densities !> @brief !!This subroutine uses the \c MPI routine \c mpi_allreduce to sum !!up the partial densities from the different nodes, using temporary !!arrays \c tmp_rho and \c tmp_current (depending on whether it !!is a scalar or vector field) in the process. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE collect_densities USE Densities, ONLY : rho,tau,current,sodens,sdens REAL(db) :: tmp_rho(nx,ny,nz,2),tmp_current(nx,ny,nz,3,2) @@ -156,15 +156,15 @@ SUBROUTINE collect_densities CALL mpi_barrier (mpi_comm_world,mpi_ierror) RETURN END SUBROUTINE collect_densities -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: Routinename !> @brief !!This subroutine collects the single-particle properties calculated !!from the wave functions and thus available only for the local wave -!!functions on each node. +!!functions on each node. !> !> @details -!!It uses a simple trick: the arrays like \c sp_energy +!!It uses a simple trick: the arrays like \c sp_energy !!are defined for the full set of indices but set to zero !!before the calculation of these properties. On each node then the !!local values are calculated but inserted at the proper index for the @@ -173,7 +173,7 @@ END SUBROUTINE collect_densities !!for each index one node contributes the correct value and the others !!zeroes. This process sounds inefficient but considering the small size !!of the arrays that does not matter. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE collect_sp_properties REAL(db) :: tmpgat(nstmax),tmpgat3(3,nstmax) CALL mpi_barrier (mpi_comm_world,mpi_ierror) @@ -203,13 +203,13 @@ SUBROUTINE collect_sp_properties sp_parity=tmpgat CALL mpi_barrier (mpi_comm_world,mpi_ierror) END SUBROUTINE collect_sp_properties -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: finish_mpi !> @brief !!This is just a wrapper for the \c MPI finalization call. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE finish_mpi - INTEGER :: ierr + INTEGER :: ierr CALL mpi_finalize(ierr) END SUBROUTINE finish_mpi END MODULE Parallel diff --git a/Code/params.f90 b/Code/params.f90 index f5dde41..97e3142 100644 --- a/Code/params.f90 +++ b/Code/params.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! MODULE: Modulename !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module contains some general parameters critical to controlling !!the code and some mathematical and physical constants used everywhere. @@ -13,8 +13,8 @@ MODULE Params !********************************************************************** !>@name General paramters !>@{ - INTEGER,PARAMETER :: db=SELECTED_REAL_KIND(12,100) !< this is a constant - !!determining the precision of real numbers in the code in a portable manner. + INTEGER,PARAMETER :: db=SELECTED_REAL_KIND(12,100) !< this is a constant + !!determining the precision of real numbers in the code in a portable manner. !!In practice it will usually be REAL(8) . ! !********************************************************************** @@ -34,15 +34,15 @@ MODULE Params !>files used by the code. Except for \c wffile, output is produced !>for an iteration or a time step at selected intervals. !>@{ - CHARACTER(LEN=80) :: wffile='none' !'NONE'. - CHARACTER(LEN=80) :: converfile='conver.res' !@} !>@name Fragment number parameters !>@{ - INTEGER,PARAMETER :: mnof=4 !@} END MODULE Params diff --git a/Code/sequential.f90 b/Code/sequential.f90 index 9071055..cd6d58b 100644 --- a/Code/sequential.f90 +++ b/Code/sequential.f90 @@ -1,17 +1,17 @@ !------------------------------------------------------------------------------ ! MODULE: Parallel !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief -!!contains the same routines as in file \c parallel.f90, but with mostly empty -!!functions to enable sequential running. +!!contains the same routines as in file \c parallel.f90, but with mostly empty +!!functions to enable sequential running. !------------------------------------------------------------------------------ MODULE Parallel USE Params, ONLY: wflag,db USE Levels, ONLY: nstmax,npsi,nstloc IMPLICIT NONE SAVE - LOGICAL,PARAMETER :: tmpi=.FALSE. + LOGICAL,PARAMETER :: tmpi=.FALSE. INTEGER, ALLOCATABLE :: node(:),localindex(:),globalindex(:) INTEGER :: mpi_nprocs,mpi_ierror,mpi_myproc, & processor_name,proc_namelen @@ -68,7 +68,7 @@ SUBROUTINE associate_nodes INTEGER :: i node=0 nstloc=nstmax - FORALL(i=1:nstmax) + FORALL(i=1:nstmax) globalindex(i)=i localindex(i)=i END FORALL diff --git a/Code/static.f90 b/Code/static.f90 index a548966..333bb2c 100644 --- a/Code/static.f90 +++ b/Code/static.f90 @@ -1,10 +1,10 @@ !------------------------------------------------------------------------------ ! MODULE: Static !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief -!!This module contains the routines needed for the static iterations. All the -!!logic needed for this case is concentrated here; all other modules except +!!This module contains the routines needed for the static iterations. All the +!!logic needed for this case is concentrated here; all other modules except !!for \c External are equally used in the dynamic calculation. !------------------------------------------------------------------------------ MODULE Static @@ -28,11 +28,11 @@ MODULE Static LOGICAL :: tvaryx_0=.FALSE. !< it .TRUE. the parameter \f$ x_0 \f$ !!is changed in every iteration in order to achieve faster convergence. INTEGER :: maxiter !< maximum number of iterations allowed. - REAL(db) :: radinx !< radius parameters (x-direction) used in the + REAL(db) :: radinx !< radius parameters (x-direction) used in the !!harmonic-oscillator initialization (see subroutine \c harmosc). - REAL(db) :: radiny !< radius parameters (y-direction) used in the + REAL(db) :: radiny !< radius parameters (y-direction) used in the !!harmonic-oscillator initialization (see subroutine \c harmosc). - REAL(db) :: radinz !< radius parameters (z-direction) used in the + REAL(db) :: radinz !< radius parameters (z-direction) used in the !!harmonic-oscillator initialization (see subroutine \c harmosc). REAL(db) :: serr !< convergence criterion. Iterations are stopped !!if the fluctuation in single-particle energies falls below this @@ -40,25 +40,25 @@ MODULE Static REAL(db) :: delesum !< used to sum up the changes in !!single-particle energies during one iteration; it is calculated in !!\c statichf but printed in \c sinfo so that it is a module variable. - REAL(db) :: x0dmp=0.2D0 !< corresponds to the parameter \f$ x_0 \f$ appearing + REAL(db) :: x0dmp=0.2D0 !< corresponds to the parameter \f$ x_0 \f$ appearing !!in the damped gradient iteration. - REAL(db) :: e0dmp=100.D0 !< corresponds to the parameter \f$ E_0 \f$ appearing + REAL(db) :: e0dmp=100.D0 !< corresponds to the parameter \f$ E_0 \f$ appearing !!in the damped gradient iteration. - REAL(db) :: x0dmpmin=0.2d0 !< corresponds to minimal value of \f$ x_0 \f$ appearing + REAL(db) :: x0dmpmin=0.2d0 !< corresponds to minimal value of \f$ x_0 \f$ appearing !!in the damped gradient iteration, if tvaryx_0 is set to .TRUE.. CONTAINS -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: getin_static !> @brief !!This routine reads the input for the static calculation using namelist -!!\c static. +!!\c static. !> !> @details !!This includes the module variables of this module, but !!also the numbers of particles, which are input quantities only in the !!static mode if initialization is not done via fragment files. In the !!case of user initialization they are also needed to correctly allocate -!!the fields; only the wave functions are then calculated in +!!the fields; only the wave functions are then calculated in !!\c user_init. The values given in the input are overwritten by !!fragment data otherwise. !! @@ -68,62 +68,62 @@ MODULE Static !!when they are computed from a formula. !! !!The variables \c charge_number and \c mass_number are also set for this case. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE getin_static NAMELIST/static/ tdiag,tlarge,maxiter, & radinx,radiny,radinz,serr,x0dmp,e0dmp,nneut,nprot,npsi,tvaryx_0 npsi=0 READ(5,static) IF(nof<=0) THEN - IF(npsi(1)==0) THEN - IF(ipair==0.OR.nof<0) THEN - npsi(1)=nneut - ELSE - npsi(1)=NINT(nneut+1.65*FLOAT(nneut)**0.666667D0) + IF(npsi(1)==0) THEN + IF(ipair==0.OR.nof<0) THEN + npsi(1)=nneut + ELSE + npsi(1)=NINT(nneut+1.65*FLOAT(nneut)**0.666667D0) IF(MOD(npsi(1),2)/=0) npsi(1)=npsi(1)+1 ENDIF ENDIF - IF(npsi(2)==0) THEN - IF(ipair==0.OR.nof<0) THEN - npsi(2)=nprot - ELSE - npsi(2)=NINT(nprot+1.65*FLOAT(nprot)**0.666667D0) + IF(npsi(2)==0) THEN + IF(ipair==0.OR.nof<0) THEN + npsi(2)=nprot + ELSE + npsi(2)=NINT(nprot+1.65*FLOAT(nprot)**0.666667D0) IF(MOD(npsi(2),2)/=0) npsi(2)=npsi(2)+1 ENDIF ENDIF - IF(nneut>npsi(1).OR.nprot>npsi(2)) & + IF(nneut>npsi(1).OR.nprot>npsi(2)) & STOP 'Particle & state numbers in conflict' nstmax=npsi(1)+npsi(2) - charge_number=nprot - mass_number=nneut+nprot + charge_number=nprot + mass_number=nneut+nprot x0dmpmin=x0dmp WRITE(*,*) "x0dmpmin=", x0dmpmin END IF END SUBROUTINE getin_static -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: init_static !> @brief !!This subroutine essentially just prints the static input and then -!!initializes the header files with their header lines. +!!initializes the header files with their header lines. !> !> @details -!!This is not included in \c getin_static, because the particle +!!This is not included in \c getin_static, because the particle !!and state numbers may have been changed by fragment input. !! -!!In addition, the damping matrices are constructed by calling +!!In addition, the damping matrices are constructed by calling !!\c setup_damping. Finally for some Skyrme forces the effective mass !!is changed to account for the center-of-mass correction. This should !!only be used if necessary and not in dynamic calculations. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE init_static IF(wflag) THEN WRITE(*,*) '***** Parameters for static calculation *****' WRITE(*,"(3(A,I4))") ' Neutrons:',nneut,' in ',npsi(1),' levels' WRITE(*,"(3(A,I4))") ' Protons :',nprot,' in ',npsi(2)-npsi(1),' levels' - WRITE(*,"(A,I6)") " Maximum number of iterations: ",maxiter + WRITE(*,"(A,I6)") " Maximum number of iterations: ",maxiter WRITE(*,"(2(A,G12.5))") ' Damping coefficient:',x0dmp, & " Damping energy scale: ",e0dmp - WRITE(*,"(A,1PE12.4)") " Convergence limit: ",serr + WRITE(*,"(A,1PE12.4)") " Convergence limit: ",serr ! initialize *.res files CALL start_protocol(converfile, & '# Iter Energy d_Energy h**2 h*h rms r^2 r^3 r^4 & @@ -143,17 +143,17 @@ SUBROUTINE init_static WRITE(*,*) '***** Nucleon mass modified for z.p.e. correction' END IF END SUBROUTINE init_static -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: statichf !> @brief !!This is the principal routine for the static iterations. It applies !!the gradient step repeatedly until convergence is achieved or the -!!maximum number of iterations is reached. +!!maximum number of iterations is reached. !> !> @details !!The following local variables !!are worth defining: -!! - \c sumflu keeps track of the fluctuations +!! - \c sumflu keeps track of the fluctuations !! in single-particle energies summed over the states. It is used as !! the convergence criterion. !! - \c addnew, \c addco are simple factors with standard @@ -178,7 +178,7 @@ END SUBROUTINE init_static !! subroutine \c skyrme is called to compute the mean-field !! components. !! - Step 3: initial gradient step: in a loop over the -!! single-particle wave functions the gradient step is applied - +!! single-particle wave functions the gradient step is applied - !! see subroutine \c grstep. The !! sum of relative changes in single-particle energies and fluctuations !! are accumulated in \c delesum and \c sumflu. This loop is @@ -204,21 +204,21 @@ END SUBROUTINE init_static !! old densities in a ratio given by \c addnew and \c addco, if !! this is turned on by \c taddnew. \c skyrme is called to !! calculate the new fields. Then the detailed single-particle -!! properties and energies are calculated and printed using +!! properties and energies are calculated and printed using !! \c sp_properties} and \c sinfo. !! - Step 9: finalizing the loop: convergence is checked by !! comparing \c sumflu per particle to \c serr, if it is smaller, !! the job terminates after writing the final wave functions. !! Otherwise, the wave functions are written if indicated by !! \c mrest. -!! . +!! . !!and the loop continues. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE statichf LOGICAL, PARAMETER :: taddnew=.TRUE. ! mix old and new densities INTEGER :: iq,nst,firstiter REAL(db) :: sumflu,denerg - REAL(db),PARAMETER :: addnew=0.2D0,addco=1.0D0-addnew + REAL(db),PARAMETER :: addnew=0.2D0,addco=1.0D0-addnew ! Step 1: initialization IF(tdiag) ALLOCATE(hmatr(nstmax,nstmax)) IF(trestart) THEN @@ -226,10 +226,10 @@ SUBROUTINE statichf ELSE iter=0 firstiter=1 - sp_energy=0.0D0 - sp_efluct1=0.0D0 + sp_energy=0.0D0 + sp_efluct1=0.0D0 sp_efluct2=0.D0 - sp_norm=0.0D0 + sp_norm=0.0D0 sumflu=0.D0 CALL schmid END IF @@ -241,19 +241,19 @@ SUBROUTINE statichf sodens=0.0D0 DO nst=1,nstmax CALL add_density(isospin(nst),wocc(nst),psi(:,:,:,:,nst), & - rho,tau,current,sdens,sodens) + rho,tau,current,sdens,sodens) ENDDO CALL skyrme ! Step 3: initial gradient step - delesum=0.0D0 - sumflu=0.0D0 + delesum=0.0D0 + sumflu=0.0D0 !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,denerg) & !$OMP SCHEDULE(STATIC) REDUCTION(+: sumflu , delesum) DO nst=1,nstmax CALL grstep(nst,isospin(nst),sp_energy(nst),denerg, & psi(:,:,:,:,nst)) - sumflu=sumflu+wocc(nst)*sp_efluct1(nst) - delesum=delesum+wocc(nst)*denerg + sumflu=sumflu+wocc(nst)*sp_efluct1(nst) + delesum=delesum+wocc(nst)*denerg ENDDO !$OMP END PARALLEL DO ! pairing and orthogonalization @@ -265,23 +265,23 @@ SUBROUTINE statichf !set x0dmp to 3* its value to get faster convergence IF(tvaryx_0) x0dmp=3.0d0*x0dmp ! step 4: start static iteration loop - Iteration: DO iter=firstiter,maxiter + Iteration: DO iter=firstiter,maxiter WRITE(*,'(a,i6,a,F12.4)') ' Static Iteration No.',iter,' x0dmp= ',x0dmp ! Step 5: gradient step - delesum=0.0D0 + delesum=0.0D0 sumflu=0.0D0 !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,denerg) & !$OMP SCHEDULE(STATIC) REDUCTION(+: sumflu , delesum) DO nst=1,nstmax CALL grstep(nst,isospin(nst),sp_energy(nst),denerg, & psi(:,:,:,:,nst)) - sumflu=sumflu+wocc(nst)*sp_efluct1(nst) - delesum=delesum+wocc(nst)*denerg + sumflu=sumflu+wocc(nst)*sp_efluct1(nst) + delesum=delesum+wocc(nst)*denerg ENDDO !$OMP END PARALLEL DO ! Step 6: diagonalize if desired IF(tdiag.AND.iter>20) THEN - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(iq) SCHEDULE(STATIC) + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(iq) SCHEDULE(STATIC) DO iq=1,2 CALL diagstep(iq,npsi(iq)-npmin(iq)+1) ENDDO @@ -304,7 +304,7 @@ SUBROUTINE statichf !$OMP REDUCTION(+:rho, tau, current, sdens, sodens) DO nst=1,nstmax CALL add_density(isospin(nst),wocc(nst),psi(:,:,:,:,nst), & - rho,tau,current,sdens,sodens) + rho,tau,current,sdens,sodens) ENDDO !$OMP END PARALLEL DO IF(taddnew) THEN @@ -323,9 +323,9 @@ SUBROUTINE statichf ! Step 9: check for convergence, saving wave functions IF(sumflu/nstmax1) THEN CALL write_wavefunctions - EXIT Iteration + EXIT Iteration END IF - IF(MOD(iter,mrest)==0) THEN + IF(MOD(iter,mrest)==0) THEN CALL write_wavefunctions ENDIF ! Step 10: update step size for the next iteration @@ -344,11 +344,11 @@ SUBROUTINE statichf END DO Iteration IF(tdiag) DEALLOCATE(hmatr) END SUBROUTINE statichf -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: grstep !> @brief !!This subroutine applies the damped gradient step to a wave -!!function \c psi. Its index is \c nst and isospin \c iq - +!!function \c psi. Its index is \c nst and isospin \c iq - !!these data are needed for the construction of the !!Hamiltonian matrix \c hmatr. The argument \c spe !!\f$ \rightarrow\epsilon \f$ represents the single-particle energy which is @@ -360,8 +360,8 @@ END SUBROUTINE statichf !! \c ps1. !! -# calculate the diagonal matrix element of \f$ \hat h \f$ !! \f[ {\tt xnormb}=\langle {\tt psin}|{\tt ps1}\rangle=\langle {\tt -!! psin}|\hat h|{\tt psin}\rangle\f] -!! and the squared norm of \c psin (in \c xnorm). The expression +!! psin}|\hat h|{\tt psin}\rangle\f] +!! and the squared norm of \c psin (in \c xnorm). The expression !! xnormb/xnorm thus corresponds !! to the expectation value \f$ \|\hat h\| \f$ in \c psin. !! @@ -387,8 +387,8 @@ END SUBROUTINE statichf !! operator acts. If \c FFT is being used for the derivatives, we use !! the routine \c laplace from module \c levels to compute !! \f[ \frac{x_{0\rm dmp}}{E_{0\rm inv}+\hat t}\,|{\tt ps1}\rangle. \f] -!! If derivatives are to done by matrices, the damping matrices -!! \c cdmpx etc. are used (see subroutine \c setdmc in module \c Grids. +!! If derivatives are to done by matrices, the damping matrices +!! \c cdmpx etc. are used (see subroutine \c setdmc in module \c Grids. !! The factors \c x0dmp and \c e0dmp have to be !! manipulated a bit in this case. !! Finally we subtract this result multiplied by \c x0act from the @@ -396,7 +396,7 @@ END SUBROUTINE statichf !! !! -# the single-particle energy is calculated from its !! new expectation value with the energy shift restore, and the -!! comparison with the initial value yields the relative change +!! comparison with the initial value yields the relative change !! \c denerg, which is passed back to the caller. !> !> @param[in] nst @@ -404,32 +404,32 @@ END SUBROUTINE statichf !> @param[in] iq !> INTEGER, takes the isospin. !> @param[in,out] spe -!> REAL(db), takes and returns single-particle energy. +!> REAL(db), takes and returns single-particle energy. !> @param[out] denerg -!> REAL(db), returns the difference in energy of the particle before and +!> REAL(db), returns the difference in energy of the particle before and !! after the iteration. !> @param[in,out] psin -!> REAL(db), array, wave function to be iterated. -!--------------------------------------------------------------------------- +!> REAL(db), array, wave function to be iterated. +!--------------------------------------------------------------------------- SUBROUTINE grstep(nst,iq,spe,denerg,psin) USE Trivial, ONLY: cmulx,cmuly,cmulz,rpsnorm,overlap !*********************************************************************** - ! - ! grstep=one damped gradient iteration step for a given - ! wave function psin with isospin iq. - ! psi=o[ psi - x0*damp*[(h-spe)psi] ] - ! + ! + ! grstep=one damped gradient iteration step for a given + ! wave function psin with isospin iq. + ! psi=o[ psi - x0*damp*[(h-spe)psi] ] + ! !*********************************************************************** INTEGER,INTENT(IN) :: nst,iq REAL(db) :: spe,denerg - COMPLEX(db) :: psin(:,:,:,:) + COMPLEX(db) :: psin(:,:,:,:) INTENT(OUT) :: denerg INTENT(INOUT) :: spe,psin REAL(db) :: x0act,esf,enrold,xnorm,xnormb,exph2,varh2 COMPLEX(db) :: ps1(nx,ny,nz,2),ps2(nx,ny,nz,2) INTEGER :: nst2 ! Step 1:(h-esf) on psin yields ps1. * - esf=spe + esf=spe CALL hpsi(iq,esf,psin,ps1) ! Step 2: calculate matrix elements xnorm=rpsnorm(psin) @@ -440,35 +440,35 @@ SUBROUTINE grstep(nst,iq,spe,denerg,psin) ENDDO IF(tdiag) hmatr(nst,nst)=hmatr(nst,nst)+spe ! Step 3: calculate fluctuation, i.e. and |h|**2 - IF((mprint>0.AND.MOD(iter,mprint)==0).OR.tvaryx_0) THEN + IF((mprint>0.AND.MOD(iter,mprint)==0).OR.tvaryx_0) THEN CALL hpsi(iq,esf,ps1,ps2) exph2=overlap(psin,ps2) varh2=rpsnorm(ps1) - sp_efluct1(nst)=SQRT(ABS(exph2/xnorm-(xnormb/xnorm)**2)) - sp_efluct2(nst)=SQRT(ABS(varh2/xnorm-(xnormb/xnorm)**2)) + sp_efluct1(nst)=SQRT(ABS(exph2/xnorm-(xnormb/xnorm)**2)) + sp_efluct2(nst)=SQRT(ABS(varh2/xnorm-(xnormb/xnorm)**2)) ENDIF ! Step 4: the damping step - IF(e0dmp>0.0D0) THEN + IF(e0dmp>0.0D0) THEN ps1=ps1 - xnormb*psin x0act=x0dmp IF(TFFT) THEN - CALL laplace(ps1,ps2,e0inv=e0dmp) + CALL laplace(ps1,ps2,e0inv=e0dmp) ELSE - CALL cmulz(cdmpz,ps1,ps2,0) - CALL cmuly(cdmpy,ps2,ps1,0) - CALL cmulx(cdmpx,ps1,ps2,0) + CALL cmulz(cdmpz,ps1,ps2,0) + CALL cmuly(cdmpy,ps2,ps1,0) + CALL cmulx(cdmpx,ps1,ps2,0) x0act=x0act/e0dmp ENDIF psin=psin - x0act*ps2 - ELSE + ELSE psin=(1.0+x0dmp*xnormb)*psin-x0dmp*ps1 ENDIF ! Step 5: energy convergence criterion - enrold=spe - spe=xnormb+esf - denerg=(enrold-spe)/ABS(spe) + enrold=spe + spe=xnormb+esf + denerg=(enrold-spe)/ABS(spe) END SUBROUTINE grstep -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: diagstep !> @brief !!This subroutine performs a diagonalization of the single-particle @@ -491,7 +491,7 @@ END SUBROUTINE grstep !! single-particle states. If \c tlarge is true, each state if !! formed independently and written onto a scratch file, otherwise an !! intermediate array \c ps1, dimensioned to hold either only -!! protons or only neutrons (number of states through the argument +!! protons or only neutrons (number of states through the argument !! \c nlin), is allocated to receive the new states, which are then !! copied back into \c psi. !> @@ -499,7 +499,7 @@ END SUBROUTINE grstep !> INTEGER, takes tisospin. !> @param[in] nlin !> INTEGER, takes the number of particles with isospin \c iq. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE diagstep(iq,nlin) !*********************************************************************** ! * @@ -568,7 +568,7 @@ END SUBROUTINE zheevd DEALLOCATE(ps1) ENDIF END SUBROUTINE diagstep -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: sinfo !> @brief !!This subroutine computes the data that are not needed for the @@ -593,26 +593,26 @@ END SUBROUTINE diagstep !!of the different contributions to the energy, a list of !!single-particle state properties, and a listing of various moments !!using \c moment_print. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE sinfo INTEGER :: il CHARACTER(*),PARAMETER :: & header=' # Par v**2 var_h1 var_h2 Norm Ekin Energy & - & Lx Ly Lz Sx Sy Sz ' + & Lx Ly Lz Sx Sy Sz ' ! calculate static observables for printout * CALL moments(0,0) CALL integ_energy CALL sum_energy ! add information to summary files - OPEN(unit=scratch,file=converfile,POSITION='APPEND') + OPEN(unit=scratch,file=converfile,POSITION='APPEND') WRITE(scratch,'(2x,i5,f9.2,3(1pg11.3),6(0pf15.8),f15.8,f15.8)') & iter,ehf,delesum/pnrtot,efluct1,efluct2,rmstot,r2tot,r3tot,r4tot,beta,gamma,avdelt(1),avdelt(2) CLOSE(scratch) - OPEN(unit=scratch,file=dipolesfile, POSITION='APPEND') + OPEN(unit=scratch,file=dipolesfile, POSITION='APPEND') WRITE(scratch,'(1x,i5,6E14.4)') iter,cmtot,cm(:,2)-cm(:,1) CLOSE(unit=scratch) - OPEN(unit=scratch,file=spinfile, POSITION='APPEND') - WRITE(scratch,'(1x,i5,9F10.4)') iter,orbital,spin,total_angmom + OPEN(unit=scratch,file=spinfile, POSITION='APPEND') + WRITE(scratch,'(1x,i5,9F10.4)') iter,orbital,spin,total_angmom CLOSE(unit=scratch) WRITE(*,'(/,A,I7,A/2(A,F12.4),A/(3(A,E12.5),A),3(A,E12.5))') & ' ***** Iteration ',iter,' *****',' Total energy: ',ehf,' MeV Total kinetic energy: ', & @@ -629,7 +629,7 @@ SUBROUTINE sinfo IF(ipair/=0) WRITE(*,'(2(A,1PE14.6))') ' Pairing energy neutrons: ', & epair(1),' protons: ',epair(2) ! output densities - IF(mplot/=0) THEN + IF(mplot/=0) THEN IF(MOD(iter,mplot)==0) THEN CALL plot_density CALL write_densities @@ -640,7 +640,7 @@ SUBROUTINE sinfo WRITE(*,'(A)') ' Neutron Single Particle States:',header DO il=1,nstmax IF(il==npmin(2)) THEN - WRITE(*,'(A)') ' Proton Single Particle States:',header + WRITE(*,'(A)') ' Proton Single Particle States:',header END IF WRITE(*,'(1X,I3,F4.0,F8.5,2F9.5,F9.6,F8.3,F10.3,3F8.3,3F7.3)') & il,sp_parity(il),wocc(il),sp_efluct1(il),sp_efluct2(il), & @@ -649,17 +649,17 @@ SUBROUTINE sinfo ENDDO CALL moment_print END SUBROUTINE sinfo -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: harmosc !> @brief -!!This subroutine initializes all wave function with harmonic oscillator states. +!!This subroutine initializes all wave function with harmonic oscillator states. !> !> @details -!!As a first step it calculates the first state, which is a Gaussian. In the -!!second step the other wave functions are obtained by multiplying the Gaussian -!!with polynomials. The subroutine does not return normalized states. In combination -!!with \c schmid the normalized lowest oscillator states are obtained. -!--------------------------------------------------------------------------- +!!As a first step it calculates the first state, which is a Gaussian. In the +!!second step the other wave functions are obtained by multiplying the Gaussian +!!with polynomials. The subroutine does not return normalized states. In combination +!!with \c schmid the normalized lowest oscillator states are obtained. +!--------------------------------------------------------------------------- SUBROUTINE harmosc USE Trivial, ONLY: rpsnorm REAL(db) :: xx,yy,zz,xx2,zz2,y2,anorm,temp @@ -670,24 +670,24 @@ SUBROUTINE harmosc !************************************************************************* ! Lowest state: Gaussian !************************************************************************* - nst=0 - DO iq=1,2 + nst=0 + DO iq=1,2 IF(iq==1) THEN nps=npsi(1) ELSE nps=npsi(2) ENDIF - ka_loop: DO ka=0,nps - DO k=0,nps - DO j=0,nps - DO i=0,nps - IF(ka==i+j+k) THEN + ka_loop: DO ka=0,nps + DO k=0,nps + DO j=0,nps + DO i=0,nps + IF(ka==i+j+k) THEN DO is=1,2 - nst=nst+1 + nst=nst+1 IF(nst>nps) EXIT ka_loop - nshell(1,nst)=i - nshell(2,nst)=j - nshell(3,nst)=k + nshell(1,nst)=i + nshell(2,nst)=j + nshell(3,nst)=k ENDDO ENDIF ENDDO @@ -698,14 +698,14 @@ SUBROUTINE harmosc END DO DO iq=1,2 nst=npmin(iq) - DO iz=1,nz + DO iz=1,nz zz2=(z(iz)/radinz)**2 - DO ix=1,nx - xx2=(x(ix)/radinx)**2 - DO iy=1,ny - y2=(y(iy)/radiny)**2 + DO ix=1,nx + xx2=(x(ix)/radinx)**2 + DO iy=1,ny + y2=(y(iy)/radiny)**2 temp=xx2+y2+zz2 - psi(ix,iy,iz,1,nst)=EXP(-(temp)) + psi(ix,iy,iz,1,nst)=EXP(-(temp)) ENDDO ENDDO ENDDO @@ -714,25 +714,25 @@ SUBROUTINE harmosc !************************************************************************* ! Higher states: lowest * polynomial !************************************************************************* - DO nst=npmin(iq)+1,npsi(iq) - is=MOD(nst-npmin(iq),2)+1 - DO iz=1,nz - IF(nshell(3,nst)/=0) THEN - zz=z(iz)**nshell(3,nst) - ELSE - zz=1.0D0 + DO nst=npmin(iq)+1,npsi(iq) + is=MOD(nst-npmin(iq),2)+1 + DO iz=1,nz + IF(nshell(3,nst)/=0) THEN + zz=z(iz)**nshell(3,nst) + ELSE + zz=1.0D0 ENDIF - DO iy=1,ny - IF(nshell(2,nst)/=0) THEN - yy=y(iy)**nshell(2,nst) - ELSE - yy=1.0D0 + DO iy=1,ny + IF(nshell(2,nst)/=0) THEN + yy=y(iy)**nshell(2,nst) + ELSE + yy=1.0D0 ENDIF - DO ix=1,nx - IF(nshell(1,nst)/=0) THEN - xx=x(ix)**nshell(1,nst) - ELSE - xx=1.0D0 + DO ix=1,nx + IF(nshell(1,nst)/=0) THEN + xx=x(ix)**nshell(1,nst) + ELSE + xx=1.0D0 ENDIF psi(ix,iy,iz,is,nst)=psi(ix,iy,iz,1,npmin(iq))*xx*yy*zz ENDDO diff --git a/Code/trivial.f90 b/Code/trivial.f90 index 825e434..97a9bed 100644 --- a/Code/trivial.f90 +++ b/Code/trivial.f90 @@ -1,7 +1,7 @@ !------------------------------------------------------------------------------ ! MODULE: Trivial !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module contains some basic calculations with wave functions and !!densities, which are similar enough to be grouped together. @@ -9,32 +9,32 @@ MODULE Trivial USE Params, ONLY: db USE Grids, ONLY: nx,ny,nz,wxyz - IMPLICIT NONE + IMPLICIT NONE CONTAINS -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: cmulx !> @brief -!!The routine does a matrix multiplication of a real matrix with a complex -!!wave function along the x-direction. +!!The routine does a matrix multiplication of a real matrix with a complex +!!wave function along the x-direction. !!\c ifadd allows accumulating results but is not used in the present code. !> !> @details !!The operation carried out here includes a loop over the spin index. !> !> @param[in] xmat -!> REAL(db), takes the matrix dimensioned +!> REAL(db), takes the matrix dimensioned !! as a square matrix with the number of points in the x-direction. !> @param[in] pinn !> COMPLEX(db), array, takes the input wave function. !> @param[out] pout !> COMPLEX(db), array, returns the multiplied wave function. !> @param[in] ifadd -!> INTEGER, if nonzero, the result of the multiplication is added to the +!> INTEGER, if nonzero, the result of the multiplication is added to the !! input wave function. -!--------------------------------------------------------------------------- - PURE SUBROUTINE cmulx(xmat,pinn,pout,ifadd) - REAL(db) :: xmat(:,:) - COMPLEX(db) :: pinn(:,:,:,:),pout(:,:,:,:) +!--------------------------------------------------------------------------- + PURE SUBROUTINE cmulx(xmat,pinn,pout,ifadd) + REAL(db) :: xmat(:,:) + COMPLEX(db) :: pinn(:,:,:,:),pout(:,:,:,:) INTEGER :: ifadd INTENT(IN) :: xmat,pinn,ifadd INTENT(INOUT) :: pout @@ -44,30 +44,30 @@ PURE SUBROUTINE cmulx(xmat,pinn,pout,ifadd) pout(ix,iy,iz,is)=pout(ix,iy,iz,is)+SUM(xmat(ix,:)*pinn(:,iy,iz,is)) END FORALL END SUBROUTINE cmulx -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: cmuly !> @brief -!!The routine does a matrix multiplication of a real matrix with a complex -!!wave function along the y-direction. +!!The routine does a matrix multiplication of a real matrix with a complex +!!wave function along the y-direction. !!\c ifadd allows accumulating results but is not used in the present code. !> !> @details !!The operation carried out here includes a loop over the spin index. !> !> @param[in] ymat -!> REAL(db), takes the matrix dimensioned +!> REAL(db), takes the matrix dimensioned !! as a square matrix with the number of points in the y-direction. !> @param[in] pinn !> COMPLEX(db), array, takes the input wave function. !> @param[out] pout !> COMPLEX(db), array, returns the multiplied wave function. !> @param[in] ifadd -!> INTEGER, if nonzero, the result of the multiplication is added to the +!> INTEGER, if nonzero, the result of the multiplication is added to the !! input wave function. -!--------------------------------------------------------------------------- - PURE SUBROUTINE cmuly(ymat,pinn,pout,ifadd) - REAL(db) :: ymat(:,:) - COMPLEX(db) :: pinn(:,:,:,:),pout(:,:,:,:) +!--------------------------------------------------------------------------- + PURE SUBROUTINE cmuly(ymat,pinn,pout,ifadd) + REAL(db) :: ymat(:,:) + COMPLEX(db) :: pinn(:,:,:,:),pout(:,:,:,:) INTEGER :: ifadd INTENT(IN) :: ymat,pinn,ifadd INTENT(INOUT) :: pout @@ -77,30 +77,30 @@ PURE SUBROUTINE cmuly(ymat,pinn,pout,ifadd) pout(ix,iy,iz,is)=pout(ix,iy,iz,is)+SUM(ymat(iy,:)*pinn(ix,:,iz,is)) END FORALL END SUBROUTINE cmuly -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: cmulz !> @brief -!!The routine does a matrix multiplication of a real matrix with a complex -!!wave function along the z-direction. +!!The routine does a matrix multiplication of a real matrix with a complex +!!wave function along the z-direction. !!\c ifadd allows accumulating results but is not used in the present code. !> !> @details !!The operation carried out here includes a loop over the spin index. !> !> @param[in] zmat -!> REAL(db), takes the matrix dimensioned +!> REAL(db), takes the matrix dimensioned !! as a square matrix with the number of points in the z-direction. !> @param[in] pinn !> COMPLEX(db), array, takes the input wave function. !> @param[out] pout !> COMPLEX(db), array, returns the multiplied wave function. !> @param[in] ifadd -!> INTEGER, if nonzero, the result of the multiplication is added to the +!> INTEGER, if nonzero, the result of the multiplication is added to the !! input wave function. -!--------------------------------------------------------------------------- - PURE SUBROUTINE cmulz(zmat,pinn,pout,ifadd) - REAL(db) :: zmat(:,:) - COMPLEX(db) :: pinn(:,:,:,:),pout(:,:,:,:) +!--------------------------------------------------------------------------- + PURE SUBROUTINE cmulz(zmat,pinn,pout,ifadd) + REAL(db) :: zmat(:,:) + COMPLEX(db) :: pinn(:,:,:,:),pout(:,:,:,:) INTEGER :: ifadd INTENT(IN) :: zmat,pinn,ifadd INTENT(INOUT) :: pout @@ -123,10 +123,10 @@ PURE SUBROUTINE cmulz(zmat,pinn,pout,ifadd) ENDDO ENDDO END SUBROUTINE cmulz -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: rpsnorm !> @brief -!!This function with one argument calculates the norm of a wave function +!!This function with one argument calculates the norm of a wave function !> !> @details !!i.e., @@ -137,13 +137,13 @@ END SUBROUTINE cmulz !> !> @param[in] ps !> COMPLEX(db), array, takes the wave function. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- PURE FUNCTION rpsnorm(ps) RESULT(r) COMPLEX(db),INTENT(IN) :: ps(:,:,:,:) REAL(db) :: r r=wxyz*SUM(CONJG(ps)*ps) END FUNCTION rpsnorm -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: overlap !> @brief !!This calculates the overlap of two wave functions \f$ psi_L \f$ and \f$ psi_R \f$ @@ -152,10 +152,10 @@ END FUNCTION rpsnorm !!This is defined as !!\f{eqnarray}{ !! \langle\psi_L|\psi_R\rangle&=&\sum_{s=\pm{1\over2}}\int \psi_L^*(\vec -!! r)\,\psi_R(\vec r)\,\D^3 r \\ +!! r)\,\psi_R(\vec r)\,\D^3 r \\ !! &\rightarrow& !! {\tt wxyz}\sum_{s=1}^2 -!! \sum_{i=1}^{nx}\sum_{j=1}^{ny}\sum_{k=1}^{nz} +!! \sum_{i=1}^{nx}\sum_{j=1}^{ny}\sum_{k=1}^{nz} !! {\tt CONJG(pl(i,j,k,s))pr(i,j,k,s)}. !!\f} !> @@ -163,117 +163,117 @@ END FUNCTION rpsnorm !> COMPEX(db), array, takes \f$ psi_L \f$. !> @param[in] pr !> COMPEX(db), array, takes \f$ psi_R \f$. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- PURE FUNCTION overlap(pl,pr) RESULT(c) COMPLEX(db) :: c,pl(:,:,:,:),pr(:,:,:,:) INTENT(IN) :: pl,pr c=wxyz*SUM(CONJG(pl)*pr) END FUNCTION overlap -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: rmulx !> @brief -!!The routine does a matrix multiplication of a real matrix with a real -!!field along the x-direction. +!!The routine does a matrix multiplication of a real matrix with a real +!!field along the x-direction. !!\c ifadd allows accumulating or substracting results. !> !> @param[in] xmat -!> REAL(db), takes the matrix dimensioned +!> REAL(db), takes the matrix dimensioned !! as a square matrix with the number of points in the x-direction. !> @param[in] finn !> REAL(db), array, takes the input wave function. !> @param[out] fout !> REAL(db), array, returns the multiplied wave function. !> @param[in] ifadd -!> INTEGER, if positive, the result of the multiplication is added to the +!> INTEGER, if positive, the result of the multiplication is added to the !! input wave function, if negative it is substracted. -!--------------------------------------------------------------------------- - PURE SUBROUTINE rmulx(xmat,finn,fout,ifadd) +!--------------------------------------------------------------------------- + PURE SUBROUTINE rmulx(xmat,finn,fout,ifadd) INTEGER :: ifadd REAL(db) :: xmat(:,:),finn(:,:,:),fout(:,:,:) INTENT(IN) :: xmat,finn,ifadd INTENT(INOUT) :: fout INTEGER :: ix,iy,iz - IF(ifadd==0) fout=0.D0 + IF(ifadd==0) fout=0.D0 IF(ifadd>=0) THEN FORALL(ix=1:nx,iy=1:ny,iz=1:nz) fout(ix,iy,iz)=fout(ix,iy,iz)+SUM(xmat(ix,:)*finn(:,iy,iz)) END FORALL - ELSE + ELSE FORALL(ix=1:nx,iy=1:ny,iz=1:nz) fout(ix,iy,iz)=fout(ix,iy,iz)-SUM(xmat(ix,:)*finn(:,iy,iz)) END FORALL ENDIF END SUBROUTINE rmulx -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: rmuly !> @brief -!!The routine does a matrix multiplication of a real matrix with a real -!!field along the y-direction. +!!The routine does a matrix multiplication of a real matrix with a real +!!field along the y-direction. !!\c ifadd allows accumulating or substracting results. !> !> @param[in] ymat -!> REAL(db), takes the matrix dimensioned +!> REAL(db), takes the matrix dimensioned !! as a square matrix with the number of points in the y-direction. !> @param[in] finn !> REAL(db), array, takes the input wave function. !> @param[out] fout !> REAL(db), array, returns the multiplied wave function. !> @param[in] ifadd -!> INTEGER, if positive, the result of the multiplication is added to the +!> INTEGER, if positive, the result of the multiplication is added to the !! input wave function, if negative it is substracted. -!--------------------------------------------------------------------------- - PURE SUBROUTINE rmuly(ymat,finn,fout,ifadd) +!--------------------------------------------------------------------------- + PURE SUBROUTINE rmuly(ymat,finn,fout,ifadd) INTEGER :: ifadd REAL(db) :: ymat(:,:),finn(:,:,:),fout(:,:,:) INTENT(IN) :: ymat,finn,ifadd INTENT(INOUT) :: fout INTEGER :: ix,iy,iz - IF(ifadd==0) fout=0.D0 - IF(ifadd>=0) THEN + IF(ifadd==0) fout=0.D0 + IF(ifadd>=0) THEN FORALL(ix=1:nx,iy=1:ny,iz=1:nz) fout(ix,iy,iz)=fout(ix,iy,iz)+SUM(ymat(iy,:)*finn(ix,:,iz)) END FORALL - ELSE + ELSE FORALL(ix=1:nx,iy=1:ny,iz=1:nz) fout(ix,iy,iz)=fout(ix,iy,iz)-SUM(ymat(iy,:)*finn(ix,:,iz)) END FORALL ENDIF END SUBROUTINE rmuly -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: rmulz !> @brief -!!The routine does a matrix multiplication of a real matrix with a real -!!field along the z-direction. +!!The routine does a matrix multiplication of a real matrix with a real +!!field along the z-direction. !!\c ifadd allows accumulating or substracting results. !> !> @param[in] zmat -!> REAL(db), takes the matrix dimensioned +!> REAL(db), takes the matrix dimensioned !! as a square matrix with the number of points in the z-direction. !> @param[in] finn !> REAL(db), array, takes the input wave function. !> @param[out] fout !> REAL(db), array, returns the multiplied wave function. !> @param[in] ifadd -!> INTEGER, if positive, the result of the multiplication is added to the +!> INTEGER, if positive, the result of the multiplication is added to the !! input wave function, if negative it is substracted. -!--------------------------------------------------------------------------- - PURE SUBROUTINE rmulz(zmat,finn,fout,ifadd) +!--------------------------------------------------------------------------- + PURE SUBROUTINE rmulz(zmat,finn,fout,ifadd) INTEGER :: ifadd REAL(db) :: zmat(:,:),finn(:,:,:),fout(:,:,:) INTENT(IN) :: zmat,finn,ifadd INTENT(INOUT) :: fout INTEGER :: ix,iy,iz - IF(ifadd==0) fout=0.D0 - IF(ifadd>=0) THEN + IF(ifadd==0) fout=0.D0 + IF(ifadd>=0) THEN FORALL(ix=1:nx,iy=1:ny,iz=1:nz) fout(ix,iy,iz)=fout(ix,iy,iz)+SUM(zmat(iz,:)*finn(ix,iy,:)) END FORALL - ELSE + ELSE FORALL(ix=1:nx,iy=1:ny,iz=1:nz) fout(ix,iy,iz)=fout(ix,iy,iz)-SUM(zmat(iz,:)*finn(ix,iy,:)) END FORALL ENDIF - RETURN + RETURN END SUBROUTINE rmulz !*********************************************************************** END MODULE Trivial diff --git a/Code/twobody.f90 b/Code/twobody.f90 index a73df53..ccb042d 100644 --- a/Code/twobody.f90 +++ b/Code/twobody.f90 @@ -1,10 +1,10 @@ !------------------------------------------------------------------------------ ! MODULE: Twobody !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief -!!This module contains the routines to distinguish between a one and two-body -!!system and to perform a two-body analysis in the case of a two-body collision +!!This module contains the routines to distinguish between a one and two-body +!!system and to perform a two-body analysis in the case of a two-body collision !!including scattering angles and momenta. !------------------------------------------------------------------------------ MODULE Twobody @@ -32,45 +32,45 @@ MODULE Twobody REAL(db) :: mass(0:2),charge(0:2),cmfrag(3,0:2),momentum(3,0:2),angmom(3,0:2) !>@} ! Data associated with the Coulomb calculation - REAL(db),ALLOCATABLE,PRIVATE :: rhotot(:,:,:) !< total density + REAL(db),ALLOCATABLE,PRIVATE :: rhotot(:,:,:) !< total density REAL(db),ALLOCATABLE,PRIVATE :: currtot(:,:,:,:) !< total current density - INTEGER,ALLOCATABLE,PRIVATE :: frag(:,:,:) !< for each mesh pint, tells whether it is in + INTEGER,ALLOCATABLE,PRIVATE :: frag(:,:,:) !< for each mesh pint, tells whether it is in !! fragment 1 or 2, or below the cutoff density, (value 0) LOGICAL :: istwobody !< principal return value. Tells whether the system is separated. CONTAINS -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: twobody_analysis !> @brief !! For a two-body collision, this subroutine checks whether the system -!! is separated and in that case determines either only the distance of +!! is separated and in that case determines either only the distance of !! the fragments or, for the final state, produces a detailed analysis !> !> @details -!! The subroutine starts off by allocating arrays \c rhotot and, if the full -!! analysis is desired, \c currtot, which are then assigned the total density and +!! The subroutine starts off by allocating arrays \c rhotot and, if the full +!! analysis is desired, \c currtot, which are then assigned the total density and !! current density, respectively. These are set to zero wherever the density !! is below the \c threshold. Then the following steps are taken: !! - For the full analysis the total current density and the isospin-dependent -!! density and current density are also cut off. Then the subroutine -!! \c analyze is called for the whole system by setting \c frag such that +!! density and current density are also cut off. Then the subroutine +!! \c analyze is called for the whole system by setting \c frag such that !! all points have fragment index 0. The coordinate values are corrected -!! for a shift in the center-of-mass. The total current density is then +!! for a shift in the center-of-mass. The total current density is then !! adjusted to remove any overall motion of the system in the grid. !! - Following this, subroutine \c cut_system is called to calculation !! fragment properties in the case of separation. If \c full_analysis is true, !! this should be the case, as it is done only for the final state. !! - For \c full_analysis=false, the analysis is finished, since the fragment !! distance has been calculated. Therefore there is a \c RETURN after the arrays -!! have been deallocated. -!! - After this the subroutine \c analyze is called separately for the two +!! have been deallocated. +!! - After this the subroutine \c analyze is called separately for the two !! fragments to determine their properties !! - Based on this information, subroutine \c scattering is called to calculate !! relative motion and the overall scattering results. !> !> @param[in] full_analysis -!> LOGICAL, if TRUE the full analysis is done, otherwise only the +!> LOGICAL, if TRUE the full analysis is done, otherwise only the !> fragment separation is calculated. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE twobody_analysis(full_analysis) LOGICAL,INTENT(IN) :: full_analysis INTEGER :: i,j,k @@ -86,7 +86,7 @@ SUBROUTINE twobody_analysis(full_analysis) FORALL(i=1:nx,j=1:ny,k=1:nz,rhotot(i,j,k) @brief !! This subroutine calculates properties of the two separate fragment !! or of the whole system. !> -!> @details +!> @details !! It evaluates integrated quantities for a subset of the mesh controlled !! by its argument \c ifrag in conjunction with the array \c !! frag. More specifically, \c ifrag=0 is used to calculate on the !! whole mesh while \c ifrag=1 or 2 select the points belonging to the -!! corresponding fragment. The array \c frag contains the relevant index +!! corresponding fragment. The array \c frag contains the relevant index !! for each mesh point. !! - Most of the calculated quantities have an index \c ifrag to store !! these results separately. Thus \c cm(1:3,0:2) is used to store the @@ -146,7 +146,7 @@ END SUBROUTINE twobody_analysis !> !> @param[in] ifrag !> INTEGER, index of fragment 1 or 2 to be investigated, or 0 for whole system -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE analyze(ifrag) INTEGER,INTENT(IN) :: ifrag REAL(db) :: velcorr(3,0:2),omega(3),inertia(3,3),curr(3),xx,yy,zz @@ -217,7 +217,7 @@ SUBROUTINE analyze(ifrag) WRITE(*,100) 'Omega from L',omega 100 FORMAT(A/(3D15.6)) END SUBROUTINE analyze -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: gauss !> @brief !! This subroutine does a simple direct solution of the 3 by 3 linear @@ -228,11 +228,11 @@ END SUBROUTINE analyze !> !> @param[in] a(3,3) !> REAL(db), Matrix of the linear system. -!> @param[in] b(3) +!> @param[in] b(3) !> REAL(db), right-hand side of the equations. !> @param[out] c(3) !> REAL(db), solution vector. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE gauss(a,b,c) REAL(db) :: a(3,3),b(3),c(3),det INTENT(IN) :: a,b @@ -253,26 +253,26 @@ SUBROUTINE gauss(a,b,c) -a(1,1)*a(2,2)*b(3))/det ENDIF END SUBROUTINE gauss - !*************************************************************** -!--------------------------------------------------------------------------- + !*************************************************************** +!--------------------------------------------------------------------------- ! DESCRIPTION: scattering !> @brief -!! This subroutine calculates the relative motion and the various energies -!! in the final state and, by referring back to the initial state before the +!! This subroutine calculates the relative motion and the various energies +!! in the final state and, by referring back to the initial state before the !! collision, also the net scattering results. !> !> @details !! The subroutine calculates the scattering angle as the sum of three contributions: -!! the ingoing Rutherford angle, the rotation angle during the TDHF calculation, +!! the ingoing Rutherford angle, the rotation angle during the TDHF calculation, !! and the outgoing Rutherford angle. -!! +!! !! The Rutherford trajectories are hyperbolas described in polar !! coordinates as \f[r(\phi)=\frac{k}{\epsilon\cos\phi-1}.\f] !! with \f$r\f$ and \f$\phi\f$ describing the relative vector of the two-body problem. !! It is confined to the angular region \f$-\phi_\infty<\phi<\phi_\infty\f$ !! with the asymptotic angle \f$\phi_\infty=\cos^{-1}(\epsilon^{-1})\f$. The !! angle \f$\phi=0\f$ corresponds to the symmetry axis of the hyperbola. -!! +!! !! A problem is that this axis is not directly connected to the !! coordinates used for the TDHF calculation, which is assumed to have !! \f$x\f$ and \f$z\f$ as the reaction plane. Since, however, the parameters of @@ -280,12 +280,12 @@ END SUBROUTINE gauss !! parameter and energy of the relative motion, it suffices to know \f$r\f$ !! to calculate the angle corresponding to it in the hyperbola system of !! coordinates. This is done slightly differently for the incoming and -!! outgoing trajectories. -!! +!! outgoing trajectories. +!! !! The following steps are done: !! - Examine the initial state before the collision: Since the properties -!! of the fragments may not be available anymore after a restart, we have to -!! partially restore them from the original input files. This requires +!! of the fragments may not be available anymore after a restart, we have to +!! partially restore them from the original input files. This requires !! that these files, or symbolic links to them, be present in current directory. !! From \c for005 it reads the center-of-mass energy \f$E_{\rm cm}\f$ and !! impact parameter \f$b\f$ as well as the initial positions of the nuclei @@ -293,10 +293,10 @@ END SUBROUTINE gauss !! \c fcent(:,i) in the code. Then the masses \f$A_i\f$ and charges \f$Z_i\f$ !! are input from the fragment data files. !! - As the next step, the variable \c mass_case is determined to make sure the -!! mass relation is the same in the ingoing and outgoing case. For symmetric +!! mass relation is the same in the ingoing and outgoing case. For symmetric !! collisions of course the scattering is not uniquely determined. !! From these data the following quantities are successively calculated: -!! -# The initial relative vector for the TDHF calculation +!! -# The initial relative vector for the TDHF calculation !! \f[\vec r(t=0)=\vec r_2(t=0)-\vec r_1(t=0),\f] !! -# the initial reduced mass !! \f[\mu=m_0\frac{A_1A_2}{A_1+A_2},\f] @@ -320,23 +320,23 @@ END SUBROUTINE gauss !! scattering angle is just the difference, !! \f[\phi_{\rm incoming}=\phi_\infty-\phi(t=0).\f] !! - The rotation of the system during the TDHF calculation is obtained from the -!! initial and final positions of the fragments via the scalar product of the +!! initial and final positions of the fragments via the scalar product of the !! relative vectors. !! - The code then calculates the relative motion properties in the final state: -!! reduced mass \c mu, relative velocity \c vrel, kinetic energy \c erel, -!! and angular momentum of relative motion \c relangmom. +!! reduced mass \c mu, relative velocity \c vrel, kinetic energy \c erel, +!! and angular momentum of relative motion \c relangmom. !! In addition the point-charge Coulomb interaction energy \c ecoul_point, !! the velocity and kinetic energy in the radial direction \c rdot and \c edot !! are produced for information. !! - For the extrapolation of the relative motion energy to infinite separation -!! a more refinde calculation of the Coulomb energy is needed. This is done by -!! using subroutine \c poisson three times: once for the total system and then -!! for each fragment individually by removing the charge density from all space points +!! a more refinde calculation of the Coulomb energy is needed. This is done by +!! using subroutine \c poisson three times: once for the total system and then +!! for each fragment individually by removing the charge density from all space points !! associated with the other fragment. This yields \c ecoul(0:2) and the Coulomb -!! interaction energy \c ecoul_int can then be obtained as +!! interaction energy \c ecoul_int can then be obtained as !! ecoul(0)-ecoul(1)-ecoul(2). The asymptotic energy is then the sum of \c erel and !! \c ecoul_int. -!! - Finally the deflection angle for the outgoing Rutherford trajectory +!! - Finally the deflection angle for the outgoing Rutherford trajectory !! is obtained. The main difference to !! the incoming trajectory is that now the quantities determining the !! hyperbola parameters are taken from the two-body analysis, namely the @@ -348,7 +348,7 @@ END SUBROUTINE gauss !! the final state. !! - The subroutine then prints all the details. !> -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE scattering ! Varaiables needed for the calculation of the initial ! Rutherford trajectory @@ -356,7 +356,7 @@ SUBROUTINE scattering REAL(db) :: fcent(3,2) ! c.m. vectors for the fragments REAL(db) :: fboost(3,2) ! boost vectors for the fragments REAL(db) :: b ! impact parameter - REAL(db) :: ecm ! relative kinetic energy + REAL(db) :: ecm ! relative kinetic energy REAL(db) :: mu ! reduced mass REAL(db) :: k,a,epsilon,p,pmax ! hyperbola parameters REAL(db) :: time ! needed for the namelist @@ -429,7 +429,7 @@ SUBROUTINE scattering roft=SQRT(SUM(final_relvec**2)) tdhf_angle=ACOS(SUM(init_relvec*final_relvec)/roft)*180.D0/pi WRITE(*,'(A,f8.3)') 'TDHF rotation angle: ',tdhf_angle - ! + ! ! Calculate relative motion quantities- ! reduced mass mu=nucleon_mass*mass(1)*mass(2)/(mass(1)+mass(2)) ! reduced mass @@ -478,7 +478,7 @@ SUBROUTINE scattering WRITE(*,'(3(A,F10.4))') 'Total Coulomb energy: ',ecoul(0), & ' Fragment 1: ',ecoul(1),' Fragment 2:',ecoul(2), & ' Coulomb interaction energy:',ecoul_int - ! + ! final_ecm=erel+ecoul_int ! now calculate final Rutherford scattering angle k=(relangmom(2)*hbc)**2/(e2*charge(1)*charge(2)*mu) @@ -504,7 +504,7 @@ SUBROUTINE scattering ' Final: ',final_ecm WRITE(*,'(3(A15,F10.3))') 'Radial: ',edot,' Centrifugal:', & e_centrifugal,' Coulomb:',ecoul_int - WRITE(*,'(2(A20,F10.3))') 'Initial L: ',initial_L,' Final orbital L: ',& + WRITE(*,'(2(A20,F10.3))') 'Initial L: ',initial_L,' Final orbital L: ',& relangmom(2) WRITE(*,'(A20,2F10.3,A20,F10.3)') 'Fragment L: ',angmom(2,1:2), & 'Final L_tot: ',relangmom(2)+angmom(2,1)+angmom(2,2) @@ -513,7 +513,7 @@ SUBROUTINE scattering WRITE(*,'(3(A15,F10.3))') 'Initial: ',initial_angle, & ' TDHF part: ',tdhf_angle,' Final: ',final_angle END SUBROUTINE scattering -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: cut_system !> @brief !! This subroutine finds out whether the system can be decomposed into two @@ -525,74 +525,74 @@ END SUBROUTINE scattering !! tensor principal axis using function \c getslope. After calculating the !! fragment centers-of-mass, this is corrected iteratively. Within the loop !! the following steps are executed: -!! -# If this is not the first iteration, calculate the slope of the +!! -# If this is not the first iteration, calculate the slope of the !! center-connecting line from the centers of mass. !! -# Call \c divpoint to find whether there is a sufficiently low density -!! between two fragments along that line. If there is, keep track of the +!! between two fragments along that line. If there is, keep track of the !! minimum-density location. -!! -# Taking the line from that point, calculate the properties of the -!! line perpendicular to the inter-center connection: the line dividing space +!! -# Taking the line from that point, calculate the properties of the +!! line perpendicular to the inter-center connection: the line dividing space !! between the fragments. !! -# Then assign fragment index \c frag(ix,iy,iz) for each grid point by -!! checking which side of the dividing line they are on. At the same time, +!! checking which side of the dividing line they are on. At the same time, !! accumulate data for the fragment masses and centers of mass. !! -# If the centers of mass have not changed significantly, stop the !! iterations. Otherwise repeat. !! -# Record the center distance and return if this is not a full analysis. !! -# Finally, make sure that the assignment of fragment index to heavier -!! or lighter fragment remains the same. This uses array \c swap to interchange +!! or lighter fragment remains the same. This uses array \c swap to interchange !! 1 and 2. !> !> @param[in] full_analysis !! LOGICAL, if false, do only what is necessary for calculating the !! distance -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- SUBROUTINE cut_system(full_analysis) LOGICAL,INTENT(IN) :: full_analysis REAL(db) :: cent(3,2),center(3,2),angle,slopev,xx,zz,diff,vol, & bb,xmin,zmin,slope INTEGER :: ix,iy,iz,itcm,ifrag INTEGER,PARAMETER :: swap(0:2)=(/ 0,2,1 /) - Iteration: DO itcm=1,10 + Iteration: DO itcm=1,10 ! calculate slope from centers of mass - IF(itcm>1) THEN - cent=center - slope=(cent(3,2)-cent(3,1))/(cent(1,2)-cent(1,1)) - bb=cent(3,1)-slope*cent(1,1) + IF(itcm>1) THEN + cent=center + slope=(cent(3,2)-cent(3,1))/(cent(1,2)-cent(1,1)) + bb=cent(3,1)-slope*cent(1,1) ELSE ! first iteration cent=0.D0 slope=getslope() bb=0.D0 ENDIF ! check whether separation - istwobody=divpoint(xmin,zmin,bb,slope) + istwobody=divpoint(xmin,zmin,bb,slope) ! determine dividing line: slope and intercept - angle=ATAN(slope) - slopev=dtan(angle+pi/2.0D0) - bb=zmin-slopev*xmin + angle=ATAN(slope) + slopev=dtan(angle+pi/2.0D0) + bb=zmin-slopev*xmin ! assign fragment cells center(:,1:2)=0.D0 - mass(1:2)=0.0D0 + mass(1:2)=0.0D0 charge(1:2)=0.0D0 - DO iz=1,nz - zz=z(iz) - DO ix=1,nx - xx=x(ix) - diff=zz-slopev*xx-bb - DO iy=1,ny - vol=wxyz + DO iz=1,nz + zz=z(iz) + DO ix=1,nx + xx=x(ix) + diff=zz-slopev*xx-bb + DO iy=1,ny + vol=wxyz frag(ix,iy,iz)=0 IF(rhotot(ix,iy,iz)mass(2)) RETURN END SELECT ! incorrect case: need to exchange indices - DO iz=1,nz - DO iy=1,ny - DO ix=1,nx + DO iz=1,nz + DO iy=1,ny + DO ix=1,nx frag(ix,iy,iz)=swap(frag(ix,iy,iz)) END DO END DO END DO END SUBROUTINE cut_system -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: getslope !> @brief -!!In this function the slope of the fragment connection line in the -!! \f$x-z \f$ plane is determined. +!!In this function the slope of the fragment connection line in the +!! \f$x-z \f$ plane is determined. !> !> @details !! The slope is determined by the @@ -662,56 +662,56 @@ END SUBROUTINE cut_system !! resulting slope is returned in the module variable \c slope. !> !> @retval newslope REAL(db), calculated slope. -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- FUNCTION getslope() RESULT(newslope) REAL(db) :: q2(3,3),xx,yy,zz,vol,denom,newslope - INTEGER :: ix,iy,iz - q2=0.D0 - DO iz=1,nz + INTEGER :: ix,iy,iz + q2=0.D0 + DO iz=1,nz zz=z(iz) - DO iy=1,ny + DO iy=1,ny yy=y(iy) - DO ix=1,nx + DO ix=1,nx xx=x(ix) vol=wxyz*rhotot(ix,iy,iz) - q2(1,1)=q2(1,1)+(xx*xx+xx*xx-yy*yy-zz*zz)*vol - q2(2,2)=q2(2,2)+(yy*yy+yy*yy-xx*xx-zz*zz)*vol - q2(3,3)=q2(3,3)+(zz*zz+zz*zz-xx*xx-yy*yy)*vol - q2(1,2)=q2(1,2)+3.D0*xx*yy*vol - q2(1,3)=q2(1,3)+3.D0*xx*zz*vol - q2(2,3)=q2(2,3)+3.D0*yy*zz*vol + q2(1,1)=q2(1,1)+(xx*xx+xx*xx-yy*yy-zz*zz)*vol + q2(2,2)=q2(2,2)+(yy*yy+yy*yy-xx*xx-zz*zz)*vol + q2(3,3)=q2(3,3)+(zz*zz+zz*zz-xx*xx-yy*yy)*vol + q2(1,2)=q2(1,2)+3.D0*xx*yy*vol + q2(1,3)=q2(1,3)+3.D0*xx*zz*vol + q2(2,3)=q2(2,3)+3.D0*yy*zz*vol ENDDO ENDDO ENDDO - q2(2,1)=q2(1,2) - q2(3,1)=q2(1,3) - q2(3,2)=q2(2,3) - denom=0.5D0*(q2(1,1)-q2(3,3)+SQRT((q2(1,1) & - -q2(3,3))**2+4.D0*q2(3,1)**2)) - IF(ABS(denom)<1.D-4) THEN - newslope=100.D0 - ELSE + q2(2,1)=q2(1,2) + q2(3,1)=q2(1,3) + q2(3,2)=q2(2,3) + denom=0.5D0*(q2(1,1)-q2(3,3)+SQRT((q2(1,1) & + -q2(3,3))**2+4.D0*q2(3,1)**2)) + IF(ABS(denom)<1.D-4) THEN + newslope=100.D0 + ELSE newslope=q2(1,3)/denom ENDIF END FUNCTION getslope -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! DESCRIPTION: divpoint !> @brief -!! The function divpoint examines the line determined by the axis of +!! The function divpoint examines the line determined by the axis of !! largest quadrupole moment to look for a suitable separation -!! point between two fragments. If one is found with sufficiently low -!! density, it returns \c TRUE and the properties of that point in its +!! point between two fragments. If one is found with sufficiently low +!! density, it returns \c TRUE and the properties of that point in its !! arguments. Otherwise it returns \c FALSE. !> !> @details !!To this end it looks at the behavior of the densities along this line. !!Since the line has no relation to the numerical grid, this is not !!trivial. -!! - First loop: Essentially it looks through the -!! \f$(x,z)\f$-plane to find points closer than half a grid spacing -!! to the desired line with equation z=slope*x+bb -!! (logical variable \c online). If the slope is larger than one, -!! i. e., if the nuclei are separating predominantly in the +!! - First loop: Essentially it looks through the +!! \f$(x,z)\f$-plane to find points closer than half a grid spacing +!! to the desired line with equation z=slope*x+bb +!! (logical variable \c online). If the slope is larger than one, +!! i. e., if the nuclei are separating predominantly in the !! \f$x\f$-direction, we need to take the equation !! x=(z-bb)/slope instead to get better resolution. To make the !! result monotonic along the line, the do loop in \f$z\f$ runs backward @@ -721,7 +721,7 @@ END FUNCTION getslope !! in arrays of length il. !! - Second loop: now the number of fragments nf !! is counted by examining this one-dimensional density curve, looking for -!! disconnected density humps above vacuum density. Where +!! disconnected density humps above vacuum density. Where !! the ``vacuum'' region starts and ends is recorded in variables !! n1 and n2. The logic is as follows: !! -# The logical variable in\_vacuum keeps track of whether @@ -743,84 +743,84 @@ END FUNCTION getslope !! \f$(x,z)\f$-plane. The code calculates the midpoint between the two !! positions and returns TRUE in this case. !> -!> @return +!> @return !> LOGICAL, indicates whether there is a fragment division or not. -!> @param[out] xmin +!> @param[out] xmin !> REAL(db), \f$x\f$-coordinate of the dividing point. -!> @param[out] zmin +!> @param[out] zmin !> REAL(db), \f$z\f$-coordinate of the dividing point. !> @param[in] bb !> REAL(db), intercept of the line joining the fragments.. !> @param[in] slope !> REAL(db), slope of the line joining the fragments. -!--------------------------------------------------------------------------- - LOGICAL FUNCTION divpoint(xmin,zmin,bb,slope) - LOGICAL :: online,two,in_vacuum +!--------------------------------------------------------------------------- + LOGICAL FUNCTION divpoint(xmin,zmin,bb,slope) + LOGICAL :: online,two,in_vacuum INTEGER :: iyy,iz1,iz2,idz,il,ix,iz,ixl(nx+ny+nz), & - izl(nx+ny+nz),i,n1,n2,nf - REAL(db) :: deltax,deltaz,xx,zz,rhol(nx+ny+nz) + izl(nx+ny+nz),i,n1,n2,nf + REAL(db) :: deltax,deltaz,xx,zz,rhol(nx+ny+nz) REAL(db),INTENT(OUT) :: xmin,zmin REAL(db),INTENT(IN) :: bb,slope REAL(db),PARAMETER :: vacuum=0.03D0 - ! - ! now the calculation of the dividing plane starts by examining the line - ! along the axis of largest quadrupole moment - ! + ! + ! now the calculation of the dividing plane starts by examining the line + ! along the axis of largest quadrupole moment + ! iyy=ny/2; deltax=0.5*dx deltaz=0.5*dz ! Loop 1 - IF(slope>=0.D0) THEN - iz1=1; iz2=nz; idz=1 - ELSE - iz1=nz; iz2=1; idz=-1 + IF(slope>=0.D0) THEN + iz1=1; iz2=nz; idz=1 + ELSE + iz1=nz; iz2=1; idz=-1 END IF - il=0 - DO ix=1,nx - xx=x(ix) - DO iz=iz1,iz2,idz - zz=z(iz) - online=.FALSE. - IF(ABS(slope)<=1.D0) THEN - online=ABS(zz-slope*xx-bb)<=deltaz - ELSE - online=ABS(xx-(zz-bb)/slope)<=deltax + il=0 + DO ix=1,nx + xx=x(ix) + DO iz=iz1,iz2,idz + zz=z(iz) + online=.FALSE. + IF(ABS(slope)<=1.D0) THEN + online=ABS(zz-slope*xx-bb)<=deltaz + ELSE + online=ABS(xx-(zz-bb)/slope)<=deltax END IF - IF(online) THEN - il=il+1 - IF(il>nx+ny+nz) THEN - WRITE(*,*) ' Increase dimensioning in function divpoint' - STOP + IF(online) THEN + il=il+1 + IF(il>nx+ny+nz) THEN + WRITE(*,*) ' Increase dimensioning in function divpoint' + STOP END IF - ixl(il)=ix; izl(il)=iz + ixl(il)=ix; izl(il)=iz rhol(il)=rhotot(ix,iyy,iz) END IF END DO END DO ! Loop 2 - nf=0; in_vacuum=.TRUE. + nf=0; in_vacuum=.TRUE. DO i=1,il - IF(rhol(i)>vacuum) THEN - IF(in_vacuum) THEN - in_vacuum=.FALSE. - nf=nf+1 - IF(nf==2) n2=i - IF(nf>2) EXIT + IF(rhol(i)>vacuum) THEN + IF(in_vacuum) THEN + in_vacuum=.FALSE. + nf=nf+1 + IF(nf==2) n2=i + IF(nf>2) EXIT END IF - ELSE - IF(.NOT.in_vacuum) THEN - IF(nf==1) n1=MAX(1,i-1) - in_vacuum=.TRUE. + ELSE + IF(.NOT.in_vacuum) THEN + IF(nf==1) n1=MAX(1,i-1) + in_vacuum=.TRUE. END IF END IF END DO - two=nf==2 + two=nf==2 ! final processing - IF(two) THEN - xmin=0.5D0*(x(ixl(n1))+x(ixl(n2))) - zmin=0.5D0*(z(izl(n1))+z(izl(n2))) + IF(two) THEN + xmin=0.5D0*(x(ixl(n1))+x(ixl(n2))) + zmin=0.5D0*(z(izl(n1))+z(izl(n2))) END IF - divpoint=two + divpoint=two END FUNCTION divpoint !*************************************************** END MODULE Twobody diff --git a/Code/user.f90 b/Code/user.f90 index 69e9053..207f483 100644 --- a/Code/user.f90 +++ b/Code/user.f90 @@ -1,11 +1,11 @@ -!--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- ! MODULE: User !> @brief !!This is included as a place to insert arbitrary user initialization of -!!the wave functions. +!!the wave functions. !> !> @details -!!It really does the same job as subroutine \c harmosc, which is a +!!It really does the same job as subroutine \c harmosc, which is a !!relatively complicated example. For this reason !!a sample user initialization is also provided in the file \c !!user_sample.f90. It produces Gaussian wave functions for three @@ -13,33 +13,33 @@ !!radii \c r. !! !!The only routine that has to be defined is \c user_init -!!\c user_init, but for more complicated initializations there can +!!\c user_init, but for more complicated initializations there can !!be any number of additional procedures accompanying it in \c user.f90. !! !!The setup assumes that the relevant particle numbers \c nneut and -!!\c nprot and numbers of states \c npmin, \c npsi, and \c nstmax are +!!\c nprot and numbers of states \c npmin, \c npsi, and \c nstmax are !!set correctly using the static input. In this case we !!assume nprot=6, nneut=6, npmin=1,7, npsi=6,12, !!and nstmax=12. Note that the occupation numbers wocc !!still have to be set explicitly, in this case they are all unity. !! -!!The routine then reads the parameters for the setup from namelist \c user. +!!The routine then reads the parameters for the setup from namelist \c user. !!This namelist can be used to read anything desired. If there !!is no user initialization, it is simply omitted from the input file. !! !!Now there is a loop over center positions with index \c ic. For !!each of them, the appropriate Gaussian is calculated and put into wave !!function # \c ic in the spin-up component; the spin-down component -!!is set to zero. Then the Gaussian is copied into index position ic+3, +!!is set to zero. Then the Gaussian is copied into index position ic+3, !!spin-down component, and finally the complete wave functions !!are copied to the proton indices by adding 6. !! -!!There is no need to orthonormalize the wave functions, since \c schmid +!!There is no need to orthonormalize the wave functions, since \c schmid !!is called before the static iterations are started. User !!initialization for the dynamic case does not appear useful; it could !!be easily done without modifying the code by running one static !!iteration and using the wave-function file generated at the beginning -!!to initialize the dynamic calculation. +!!to initialize the dynamic calculation. !--------------------------------------------------------------------------- MODULE User USE Params diff --git a/Code/ylm.f90 b/Code/ylm.f90 index 82b5de3..4cd86f2 100644 --- a/Code/ylm.f90 +++ b/Code/ylm.f90 @@ -1,11 +1,11 @@ !------------------------------------------------------------------------------ ! MODULE: Spherical Harmonics !------------------------------------------------------------------------------ -! DESCRIPTION: +! DESCRIPTION: !> @brief !!This module contains the definition of the spherical harmonics function \f$Y_{LM}\f$ !> -!>@details +!>@details !! It contains the definition \f$Y_{LM}\f$ in terms of Associate Legendre Polinomial (\f$P_{LM}\f$) Module Spherical_Harmonics USE Params, ONLY: db,pi @@ -13,8 +13,8 @@ Module Spherical_Harmonics ! DESCRIPTION: Fact(n) !> @brief !!Function Fact Calculates the Factorial. - REAL(db) function Fact(n)result(Fct) - + REAL(db) function Fact(n)result(Fct) + IMPLICIT NONE integer:: n,i ! real(db)::Fct @@ -26,8 +26,8 @@ REAL(db) function Fact(n)result(Fct) ! DESCRIPTION: Fact2(n) !> @brief !!Function Fact2 Calculates the double Factorial. - REAL(db) function Fact2(n)result(Fct2) - + REAL(db) function Fact2(n)result(Fct2) + IMPLICIT NONE integer :: n,i ! real(db)::Fct2 @@ -41,7 +41,7 @@ REAL(db) function Fact2(n)result(Fct2) ! DESCRIPTION: Plm !> @brief !!Function Plm Calculates the Associate Legendre Polinomial - real(db) function Plm(x,l,m)result(P_lm) + real(db) function Plm(x,l,m)result(P_lm) IMPLICIT NONE integer, intent(in) :: l,m integer :: em_i,em_i_lower,em_i_lower2,el,em @@ -57,7 +57,7 @@ real(db) function Plm(x,l,m)result(P_lm) P00=1 P10=x P11 = -sqrt(1-x**2) - + ALLOCATE(Ps(l+1,(2*l+1))) Ps=0.0D0 Ps(1,1)=P00 @@ -65,7 +65,7 @@ real(db) function Plm(x,l,m)result(P_lm) Ps(2,2) = P10 Ps(2,3) = -Ps(2,1)/2 do el = 2,l - Ps(el+1,el-el+1) = ((-1)**el)*Fact2(2*el-1)*(1-x**2)**(el/2) + Ps(el+1,el-el+1) = ((-1)**el)*Fact2(2*el-1)*(1-x**2)**(el/2) do em = el-1,-el,-1 em_i = el-em+1 if (em>=0)then @@ -82,7 +82,7 @@ real(db) function Plm(x,l,m)result(P_lm) ! DESCRIPTION: Y_lm !> @brief !!Function Y_lm Calculates the spherical harmonics (Real part) - real(db) function Y_lm(l,m,x,y,z)result(Ylm) + real(db) function Y_lm(l,m,x,y,z)result(Ylm) IMPLICIT NONE integer, intent(in) :: l,m real(db), intent(in) :: x,y,z @@ -107,7 +107,7 @@ real(db) function Y_lm(l,m,x,y,z)result(Ylm) Ylm = ylm_re ! write(*,*) 'positive m',l,m,abs(m),mphi*180/pi,acos(cos_theta)*180/pi,x,y,z,cos(mphi) end if - + end function diff --git a/Test/GR/spectral_analysis.f90 b/Test/GR/spectral_analysis.f90 index aa17cb7..f766bb3 100644 --- a/Test/GR/spectral_analysis.f90 +++ b/Test/GR/spectral_analysis.f90 @@ -1,15 +1,15 @@ PROGRAM spectral_analysis ! -! This program computes the spectral distribution of quadrupole -! oscillations produced as output on file 'quadrupoles.res'. +! This program computes the spectral distribution of quadrupole +! oscillations produced as output on file 'quadrupoles.res'. ! Isoscalar and isovector quadrupole moments are composed from the ! information on x**2, y**2, and z**2 in 'quadrupoles.res'. ! A limiting time profile of type COS**n is multiplied to the signal ! to guarantee that the signal approaches zero at the end of the ! time interval (called time filtering, or windowing). ! It is assumed that the quadrupole motion was initialized by an -! instantaneous isoscalar quadrupole boost. The spectral strength is -! then simply the imaginary part of the Fourier transform of the +! instantaneous isoscalar quadrupole boost. The spectral strength is +! then simply the imaginary part of the Fourier transform of the ! isoscalar signal. The code also produces the Fourier transform ! of the isovector signal. This has to be taken with care because ! the excitation operator and the analzing operator are not the same @@ -31,7 +31,7 @@ PROGRAM spectral_analysis ! -! A first quick reading to check the length and integrity of +! A first quick reading to check the length and integrity of ! the data file. ! OPEN(1,file='extfield.res') @@ -62,7 +62,7 @@ PROGRAM spectral_analysis READ(1,'(1x)') DO nt=1,ndim READ(1,*,ERR=29,END=29) timax,qext(nt) - filter = COS((nt-1)*PI/(2D0*(ndim-1)))**nfilter + filter = COS((nt-1)*PI/(2D0*(ndim-1)))**nfilter qfilter(nt) = qext(nt)*filter WRITE(2,'(f10.2,4(1pg13.5))') (nt-1)*deltime,qfilter(nt) END DO diff --git a/Utils/Cuts/Cuts.f90 b/Utils/Cuts/Cuts.f90 index fce6620..d350135 100755 --- a/Utils/Cuts/Cuts.f90 +++ b/Utils/Cuts/Cuts.f90 @@ -49,7 +49,7 @@ SUBROUTINE writeit(a,b,c,n,m,name) WRITE(*,*) 'Writing ',filename OPEN(UNIT=11,FILE=filename,FORM='FORMATTED',STATUS='REPLACE') DO j=1,m - DO i=1,n + DO i=1,n WRITE(11,'(2F8.3,E14.6)') b(i),c(j),a(i,j) END DO WRITE(11,'(1X)') diff --git a/Utils/Fileinfo/fileinfo.f90 b/Utils/Fileinfo/fileinfo.f90 index 719d49d..f031a30 100644 --- a/Utils/Fileinfo/fileinfo.f90 +++ b/Utils/Fileinfo/fileinfo.f90 @@ -12,7 +12,7 @@ PROGRAM Fileinfo OPEN(UNIT=10,FILE=filename,FORM='UNFORMATTED') ! Try to read as wave function file first READ(10,ERR=3,END=3) iter,time,forcename,nstmax,nneut,nprot,number,npsi, & - charge_number,mass_number,cm + charge_number,mass_number,cm READ(10) nx,ny,nz,dx,dy,dz,wxyz WRITE(*,'(3A,I9,A,F10.2)') 'File ',filename,'is a wave function file at iteration', & iter,' or time ',time diff --git a/Utils/Overlap/overlap.f90 b/Utils/Overlap/overlap.f90 index e424bd8..2a1a313 100755 --- a/Utils/Overlap/overlap.f90 +++ b/Utils/Overlap/overlap.f90 @@ -20,7 +20,7 @@ PROGRAM Analyze OPEN(UNIT=11,FILE=file1,STATUS ='old',FORM='unformatted') OPEN(UNIT=12,FILE=file2,STATUS ='old',FORM='unformatted') READ(11) iter,time,forcename,nstmax,nneut,nprot,number,npsi, & - charge_number,mass_number,cm + charge_number,mass_number,cm READ(12) iter,time,forcename2,nstmax2,nneut2,nprot2,number,npsi, & charge_number,mass_number,cm2 IF(nstmax/=nstmax2.OR.nneut/=nneut2.OR.nprot/=nprot2) THEN From 5a822408fb849f68ea36ef3c52565f7ed1b50c47 Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Wed, 17 Sep 2025 15:20:00 +0100 Subject: [PATCH 03/17] Fix Fortitude rule MOD011 using automatic fix flag --- Code/external.f90 | 4 ++-- Code/fourier.f90 | 12 ++++++------ Code/grids.f90 | 2 +- Code/inout.f90 | 4 ++-- Code/moment.f90 | 2 +- Code/twobody.f90 | 2 +- Code/user_sample.f90 | 2 +- Utils/Convert/Tdhf2Silo.f90 | 2 +- 8 files changed, 15 insertions(+), 15 deletions(-) diff --git a/Code/external.f90 b/Code/external.f90 index 17a1f9e..2e12920 100644 --- a/Code/external.f90 +++ b/Code/external.f90 @@ -113,8 +113,8 @@ END SUBROUTINE getin_external SUBROUTINE init_external REAL(db) :: facn,facp,facr,xlim,ylim,zlim,dip_f INTEGER :: ix,iy,iz - CHARACTER(14),PARAMETER :: pulsetype(0:2)=(/ 'Instantaneous ', & - 'Gaussian ','Cosine squared' /) + CHARACTER(14),PARAMETER :: pulsetype(0:2)=[ 'Instantaneous ', & + 'Gaussian ','Cosine squared' ] IF(ipulse<0.OR.ipulse>2) STOP & ' External field: called with invalid pulse type' IF(wflag) THEN diff --git a/Code/fourier.f90 b/Code/fourier.f90 index fcffb7c..6f57340 100644 --- a/Code/fourier.f90 +++ b/Code/fourier.f90 @@ -83,22 +83,22 @@ SUBROUTINE init_fft FFTW_FORWARD, FFTW_planflag) CALL dfftw_plan_dft_3d(pbackward,nx,ny,nz,p(:,:,:,1,1),p(:,:,:,1,1), & FFTW_BACKWARD, FFTW_planflag) - CALL dfftw_plan_many_dft(xforward,1,(/nx/),2*ny*nz, & + CALL dfftw_plan_many_dft(xforward,1,[nx],2*ny*nz, & p(:,:,:,:,1),0,1,nx,p(:,:,:,:,2),0,1,nx, & FFTW_FORWARD, FFTW_planflag) - CALL dfftw_plan_many_dft(xbackward,1,(/nx/),2*ny*nz, & + CALL dfftw_plan_many_dft(xbackward,1,[nx],2*ny*nz, & p(:,:,:,:,1),0,1,nx,p(:,:,:,:,1),0,1,nx, & FFTW_BACKWARD, FFTW_planflag) - CALL dfftw_plan_many_dft(yforward,1,(/ny/),nx, & + CALL dfftw_plan_many_dft(yforward,1,[ny],nx, & p(:,:,:,:,1),0,nx,1,p(:,:,:,:,2),0,nx,1, & FFTW_FORWARD, FFTW_planflag) - CALL dfftw_plan_many_dft(ybackward,1,(/ny/),nx, & + CALL dfftw_plan_many_dft(ybackward,1,[ny],nx, & p(:,:,:,:,1),0,nx,1,p(:,:,:,:,1),0,nx,1, & FFTW_BACKWARD, FFTW_planflag) - CALL dfftw_plan_many_dft(zforward,1,(/nz/),nx*ny, & + CALL dfftw_plan_many_dft(zforward,1,[nz],nx*ny, & p(:,:,:,:,1),0,nx*ny,1,p(:,:,:,:,2),0,nx*ny,1, & FFTW_FORWARD, FFTW_planflag) - CALL dfftw_plan_many_dft(zbackward,1,(/nz/),nx*ny, & + CALL dfftw_plan_many_dft(zbackward,1,[nz],nx*ny, & p(:,:,:,:,1),0,nx*ny,1,p(:,:,:,:,1),0,nx*ny,1, & FFTW_BACKWARD, FFTW_planflag) IF(wflag) WRITE(*,*) '***** FFTW3 plans established *****' diff --git a/Code/grids.f90 b/Code/grids.f90 index db3c589..07f4ae7 100644 --- a/Code/grids.f90 +++ b/Code/grids.f90 @@ -144,7 +144,7 @@ SUBROUTINE init_coord(name,nv,dv,v,der1v,der2v,cdmpv) der2v(:,:),cdmpv(:,:) INTEGER :: i ALLOCATE(v(nv),der1v(nv,nv),der2v(nv,nv),cdmpv(nv,nv)) - v=(/ ((i-1)*dv-0.5D0*FLOAT(nv-1)*dv,i=1,nv) /) + v=[ ((i-1)*dv-0.5D0*FLOAT(nv-1)*dv,i=1,nv) ] IF(wflag) THEN WRITE(*,'(1X,A,I3,A,F8.4,2(A,F8.4))') name // ' direction: ',nv, & ' points, spacing:',dv,' ranging from ',v(1),' to ',v(nv) diff --git a/Code/inout.f90 b/Code/inout.f90 index 6e3d09c..79a9fba 100644 --- a/Code/inout.f90 +++ b/Code/inout.f90 @@ -273,9 +273,9 @@ SUBROUTINE plot_density INTEGER,PARAMETER :: ixsc=10,izsc=6 REAL(db) :: rhoplt(nx,nz),xperi CHARACTER(1) :: ifun(121),ibor(121) - CHARACTER(1),PARAMETER :: icar(25)=(/' ','1',' ','3',' ', & + CHARACTER(1),PARAMETER :: icar(25)=[' ','1',' ','3',' ', & '5',' ','7',' ','9',' ','b',' ','d',' ','f',' ','h', & - ' ','j',' ','l',' ','n','*' /) + ' ','j',' ','l',' ','n','*' ] REAL(db) :: xco(nx+1),zco(nz+1),dimx,dimz,dxp,dzp, & xcu,zcu INTEGER :: nhor,nver,ntkx,ntkz,jcar,ntz,i,j diff --git a/Code/moment.f90 b/Code/moment.f90 index e148686..f2b7e35 100644 --- a/Code/moment.f90 +++ b/Code/moment.f90 @@ -316,7 +316,7 @@ END SUBROUTINE moment_shortprint !--------------------------------------------------------------------------- SUBROUTINE moment_print INTEGER :: iq - CHARACTER(11),PARAMETER :: Name(2)=(/ ' Neutron: ',' Proton: '/) + CHARACTER(11),PARAMETER :: Name(2)=[ ' Neutron: ',' Proton: '] Write(*,'(A)') ' Part.Num. rms-radius q20 & & & & ' diff --git a/Code/twobody.f90 b/Code/twobody.f90 index ccb042d..c6abcbe 100644 --- a/Code/twobody.f90 +++ b/Code/twobody.f90 @@ -552,7 +552,7 @@ SUBROUTINE cut_system(full_analysis) REAL(db) :: cent(3,2),center(3,2),angle,slopev,xx,zz,diff,vol, & bb,xmin,zmin,slope INTEGER :: ix,iy,iz,itcm,ifrag - INTEGER,PARAMETER :: swap(0:2)=(/ 0,2,1 /) + INTEGER,PARAMETER :: swap(0:2)=[ 0,2,1 ] Iteration: DO itcm=1,10 ! calculate slope from centers of mass IF(itcm>1) THEN diff --git a/Code/user_sample.f90 b/Code/user_sample.f90 index 8b6debd..6cf8360 100644 --- a/Code/user_sample.f90 +++ b/Code/user_sample.f90 @@ -9,7 +9,7 @@ SUBROUTINE init_user REAL(db) :: center(3),d,r NAMELIST /user/ d,r READ(5,user) - center=(/-d,0.D0,d/) + center=[-d,0.D0,d] WRITE(*,'(A,F5.2,A,F5.4)') '3 alpha cluster initialization with distance ' & ,d,' and radii ',r wocc=1.D0 diff --git a/Utils/Convert/Tdhf2Silo.f90 b/Utils/Convert/Tdhf2Silo.f90 index d11add7..7087c33 100755 --- a/Utils/Convert/Tdhf2Silo.f90 +++ b/Utils/Convert/Tdhf2Silo.f90 @@ -33,7 +33,7 @@ PROGRAM Tdhf2Silo ! put this into silo file ! ! set dimensions for following zone-centered arrays - dims=(/nx,ny,nz/) + dims=[nx,ny,nz] ! make options to add cycle number and time ret=dbmkoptlist(5,op) CALL errchk(ret,' could not make optlist') From bbf3aa22110dfe133ba40b4886da1ad3ae4cd2c0 Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Wed, 17 Sep 2025 15:20:34 +0100 Subject: [PATCH 04/17] Fix Fortitude rule MOD021 using automatic fix flag --- Code/external.f90 | 10 +++++----- Code/moment.f90 | 2 +- Code/ylm.f90 | 4 ++-- Utils/Overlap/overlap.f90 | 2 +- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/Code/external.f90 b/Code/external.f90 index 2e12920..a5c144d 100644 --- a/Code/external.f90 +++ b/Code/external.f90 @@ -140,7 +140,7 @@ SUBROUTINE init_external facn=1.0D0 facp=1.0D0 ELSE - if (L_val .eq. 1)then + if (L_val == 1)then facn=-(charge_number/mass_number) facp=(mass_number-charge_number)/mass_number else @@ -158,13 +158,13 @@ SUBROUTINE init_external DO ix=1,nx ! damped version facr = ampl_ext*SQRT(2.0d0*L_val+1.0d0) - if (L_val .ge. 2)then + if (L_val >= 2)then facr=facr*Y_lm(L_val,M_val,x(ix),y(iy),z(iz)) facr=facr*SQRT(x(ix)**2+y(iy)**2+z(iz)**2)**(L_val) - else if (L_val .eq. 0)then + else if (L_val == 0)then facr=facr*(0.5d0*SQRT(1.0d0/PI)) facr=facr*SQRT(x(ix)**2+y(iy)**2+z(iz)**2)**(2) - else if (L_val .eq. 1)then + else if (L_val == 1)then ! vol=wxyz*rho(ix,iy,iz,iq) ! write(*,*)'RMS value of r for dipole case',r_avg IF(isoext==0) THEN @@ -181,7 +181,7 @@ SUBROUTINE init_external end if facr=facr/(1.0D0+EXP((SQRT(x(ix)**2+y(iy)**2+z(iz)**2)-radext)/widext)) !< Damping is done using parameters radext and widext - if (only_P.eq.1)then + if (only_P==1)then extfield(ix,iy,iz,1)=0.0d0 extfield(ix,iy,iz,2)=facr*facp else diff --git a/Code/moment.f90 b/Code/moment.f90 index f2b7e35..3ad6e1b 100644 --- a/Code/moment.f90 +++ b/Code/moment.f90 @@ -233,7 +233,7 @@ SUBROUTINE moments(L_val,M_val) facp=1.0D0 Di_tot(iii) = facn*Di_is(1)+facp*Di_is(2) ELSE - if (L_val .eq. 1)then + if (L_val == 1)then facn=-(pnr(2)/pnrtot) facp=(pnr(1))/pnrtot Di_tot(iii) = facn*Di_iv(1)+facp*Di_iv(2) diff --git a/Code/ylm.f90 b/Code/ylm.f90 index 4cd86f2..762e84b 100644 --- a/Code/ylm.f90 +++ b/Code/ylm.f90 @@ -91,12 +91,12 @@ real(db) function Y_lm(l,m,x,y,z)result(Ylm) r=sqrt(x**2+y**2+z**2) cos_theta = z/r - if (y .ge. 0)then + if (y >= 0)then mphi = abs(m)*acos(x/(sqrt(x**2+y**2))) else mphi = -1*abs(m)*acos(x/(sqrt(x**2+y**2))) end if - if (m .lt. 0) then + if (m < 0) then const = sqrt(((2*l+1)*(Fact(l-abs(m))))/(4*PI*Fact(l+abs(m)))) ylm_im = const*Plm(cos_theta,l,abs(m))*sin(mphi) Ylm = ylm_im!*((-1)**m) diff --git a/Utils/Overlap/overlap.f90 b/Utils/Overlap/overlap.f90 index 2a1a313..a85f66b 100755 --- a/Utils/Overlap/overlap.f90 +++ b/Utils/Overlap/overlap.f90 @@ -116,7 +116,7 @@ COMPLEX(8) function determinant(a,n) COMPLEX(8) a(n,n),det(2),work(n) integer ipvt(n),info call zgefa(a,n,n,ipvt,info) - if(info.ne.0) then + if(info/=0) then write(*,*) ' Determinant error ',info stop endif From 705110ffb5e8bf010b46c64625c541d50b84acaf Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Wed, 17 Sep 2025 15:20:57 +0100 Subject: [PATCH 05/17] Fix Fortitude rule S071 using automatic fix flag --- Code/inout.f90 | 4 ++-- Utils/Convert/Tdhf2Silo.f90 | 4 ++-- Utils/Overlap/overlap.f90 | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Code/inout.f90 b/Code/inout.f90 index 79a9fba..6679fb7 100644 --- a/Code/inout.f90 +++ b/Code/inout.f90 @@ -220,7 +220,7 @@ SUBROUTINE write_one_density(name,values) CHARACTER(*),INTENT(IN) :: name REAL(db),INTENT(IN) :: values(nx,ny,nz,2) CHARACTER(10) :: stored_name - REAL(db) a(nx,ny,nz) + REAL(db) :: a(nx,ny,nz) stored_name=name WRITE(scratch) stored_name,.FALSE.,write_isospin IF(write_isospin) THEN @@ -245,7 +245,7 @@ SUBROUTINE write_vec_density(name,values) CHARACTER(*),INTENT(IN) :: name REAL(db),INTENT(IN) :: values(nx,ny,nz,3,2) CHARACTER(10) :: stored_name - REAL(db) a(nx,ny,nz,3) + REAL(db) :: a(nx,ny,nz,3) stored_name=name WRITE(scratch) stored_name,.TRUE.,write_isospin IF(write_isospin) THEN diff --git a/Utils/Convert/Tdhf2Silo.f90 b/Utils/Convert/Tdhf2Silo.f90 index 7087c33..05916ba 100755 --- a/Utils/Convert/Tdhf2Silo.f90 +++ b/Utils/Convert/Tdhf2Silo.f90 @@ -11,8 +11,8 @@ PROGRAM Tdhf2Silo CHARACTER(11) :: filename CHARACTER :: stored_name*10 INTEGER :: is,ret,iret,dims(3),dbid,op,i,ndefs - INTEGER err, ierr, types(maxdefs), lnames(maxdefs), ldefs(maxdefs) - INTEGER oldlen + INTEGER :: err, ierr, types(maxdefs), lnames(maxdefs), ldefs(maxdefs) + INTEGER :: oldlen ! Initialize some 20 character length strings CHARACTER(40) :: names(maxdefs),defs(maxdefs) ! Store the length of each string diff --git a/Utils/Overlap/overlap.f90 b/Utils/Overlap/overlap.f90 index a85f66b..e02c765 100755 --- a/Utils/Overlap/overlap.f90 +++ b/Utils/Overlap/overlap.f90 @@ -113,8 +113,8 @@ COMPLEX(8) FUNCTION scalprod(a,b) END FUNCTION scalprod COMPLEX(8) function determinant(a,n) INTEGER,INTENT(IN) :: n - COMPLEX(8) a(n,n),det(2),work(n) - integer ipvt(n),info + COMPLEX(8) :: a(n,n),det(2),work(n) + integer :: ipvt(n),info call zgefa(a,n,n,ipvt,info) if(info/=0) then write(*,*) ' Determinant error ',info From 364edf1ba7611e3e43aaa7486f5b7e4c9ecb040a Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Wed, 17 Sep 2025 15:21:23 +0100 Subject: [PATCH 06/17] Fix Fortitude rule S061 using automatic fix flag --- Code/ylm.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Code/ylm.f90 b/Code/ylm.f90 index 762e84b..17e6e37 100644 --- a/Code/ylm.f90 +++ b/Code/ylm.f90 @@ -22,7 +22,7 @@ REAL(db) function Fact(n)result(Fct) do i =1,n Fct = Fct*i end do - end function + end function Fact ! DESCRIPTION: Fact2(n) !> @brief !!Function Fact2 Calculates the double Factorial. @@ -36,7 +36,7 @@ REAL(db) function Fact2(n)result(Fct2) Fct2 = Fct2*n n=n-2 end do - end function + end function Fact2 ! DESCRIPTION: Plm !> @brief @@ -78,7 +78,7 @@ real(db) function Plm(x,l,m)result(P_lm) end do end do P_lm = Ps(l+1,l+1-m) - end function + end function Plm ! DESCRIPTION: Y_lm !> @brief !!Function Y_lm Calculates the spherical harmonics (Real part) @@ -109,6 +109,6 @@ real(db) function Y_lm(l,m,x,y,z)result(Ylm) end if - end function + end function Y_lm End Module Spherical_Harmonics \ No newline at end of file From b57533a197c62bae5b3a08e09a76df45e3f8e0d6 Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Wed, 17 Sep 2025 15:21:42 +0100 Subject: [PATCH 07/17] Fix Fortitude rule OB061 using automatic fix flag --- Utils/Convert/Tdhf2Silo.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Utils/Convert/Tdhf2Silo.f90 b/Utils/Convert/Tdhf2Silo.f90 index 05916ba..6409934 100755 --- a/Utils/Convert/Tdhf2Silo.f90 +++ b/Utils/Convert/Tdhf2Silo.f90 @@ -107,7 +107,7 @@ PROGRAM Tdhf2Silo !******************************************************************* CONTAINS SUBROUTINE write_density(rname,data) - CHARACTER*(*),INTENT(IN) :: rname + CHARACTER(len=*),INTENT(IN) :: rname REAL(8),INTENT(IN) :: data(*) WRITE(*,*) 'Writing ',rname,' ',LEN(rname) ret=dbputqv1(dbid,rname,LEN(rname),'Mesh',4,data,dims,3, & From f8f0336c4920d0470169947a10e16bc91348c15a Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Thu, 18 Sep 2025 16:02:37 +0100 Subject: [PATCH 08/17] Fix Fortitude rule S001 --- Code/external.f90 | 3 ++- Code/forces.f90 | 3 ++- Code/moment.f90 | 4 ++-- Code/static.f90 | 3 ++- Code/ylm.f90 | 5 +++-- 5 files changed, 11 insertions(+), 7 deletions(-) diff --git a/Code/external.f90 b/Code/external.f90 index a5c144d..ca254e2 100644 --- a/Code/external.f90 +++ b/Code/external.f90 @@ -171,7 +171,8 @@ SUBROUTINE init_external facr=facr*Y_lm(L_val,M_val,x(ix),y(iy),z(iz)) dip_f = (r_avg**2)*5.0d0/3.0d0 ! write(*,*)'__DIPOLE FACTOR__',r_avg,dip_f - facr=facr*(SQRT(x(ix)**2+y(iy)**2+z(iz)**2)**3 - dip_f*SQRT(x(ix)**2+y(iy)**2+z(iz)**2)) + facr=facr*(SQRT(x(ix)**2+y(iy)**2+z(iz)**2)**3 - & + dip_f*SQRT(x(ix)**2+y(iy)**2+z(iz)**2)) else if (isoext==1) then facr=facr*Y_lm(L_val,M_val,x(ix),y(iy),z(iz)) diff --git a/Code/forces.f90 b/Code/forces.f90 index 13b0983..0b6d902 100644 --- a/Code/forces.f90 +++ b/Code/forces.f90 @@ -211,7 +211,8 @@ SUBROUTINE read_force WRITE(*,"(5(A6,F12.5))") "x0",f%x0,"x1",f%x1,"x2",f%x2,"x3",f%x3,"b4p",f%b4p WRITE(*,"(A6,F12.5)") "Power",f%power WRITE(*,"(A,I2)") " Pairing parameters: Option ipair:",ipair - WRITE(*,"(3(A7,F12.5))") "v0prot",p%v0prot,"v0neut",p%v0neut,"rho0pr",p%rho0pr,'MIX',p%mixture + WRITE(*,"(3(A7,F12.5))") "v0prot",p%v0prot,"v0neut",p%v0neut,"rho0pr",p%rho0pr,'MIX',& + p%mixture ENDIF END SUBROUTINE read_force END MODULE Forces diff --git a/Code/moment.f90 b/Code/moment.f90 index 3ad6e1b..f7da666 100644 --- a/Code/moment.f90 +++ b/Code/moment.f90 @@ -193,8 +193,8 @@ SUBROUTINE moments(L_val,M_val) eta = vol*(SUM(x2))*5.0d0/3.0d0 - Di_is(iq)=vol*(Y_lm(1,M_val,x(ix),y(iy),z(iz))*(SQRT(x(ix)**2+y(iy)**2+z(iz)**2)**3- & - SQRT(x(ix)**2+y(iy)**2+z(iz)**2)*eta)*SQRT(2*1+1.0d0))+Di_is(iq) + Di_is(iq)=vol*(Y_lm(1,M_val,x(ix),y(iy),z(iz))*(SQRT(x(ix)**2+y(iy)**2+z(iz)**2)**3& + - SQRT(x(ix)**2+y(iy)**2+z(iz)**2)*eta)*SQRT(2*1+1.0d0))+Di_is(iq) ! print*,'inside loop 1' tmp=(Y_lm(1,M_val,x(ix),y(iy),z(iz))*SQRT(3.0d0)) diff --git a/Code/static.f90 b/Code/static.f90 index 333bb2c..56e60fa 100644 --- a/Code/static.f90 +++ b/Code/static.f90 @@ -606,7 +606,8 @@ SUBROUTINE sinfo ! add information to summary files OPEN(unit=scratch,file=converfile,POSITION='APPEND') WRITE(scratch,'(2x,i5,f9.2,3(1pg11.3),6(0pf15.8),f15.8,f15.8)') & - iter,ehf,delesum/pnrtot,efluct1,efluct2,rmstot,r2tot,r3tot,r4tot,beta,gamma,avdelt(1),avdelt(2) + iter,ehf,delesum/pnrtot,efluct1,efluct2,rmstot,r2tot,r3tot,r4tot,beta,gamma,avdelt(1),& + avdelt(2) CLOSE(scratch) OPEN(unit=scratch,file=dipolesfile, POSITION='APPEND') WRITE(scratch,'(1x,i5,6E14.4)') iter,cmtot,cm(:,2)-cm(:,1) diff --git a/Code/ylm.f90 b/Code/ylm.f90 index 17e6e37..c272454 100644 --- a/Code/ylm.f90 +++ b/Code/ylm.f90 @@ -71,7 +71,8 @@ real(db) function Plm(x,l,m)result(P_lm) if (em>=0)then em_i_lower = el-em em_i_lower2 = el-1-em - Ps(el+1,em_i) = (x*(2*el-1)*Ps(el,em_i_lower)-(el+em-1)*Ps(el-1,em_i_lower2))/(el-em) + Ps(el+1,em_i) = (x*(2*el-1)*Ps(el,em_i_lower)-(el+em-1)*Ps(el-1,em_i_lower2))/& + (el-em) else if (em<0)then Ps(el+1, em_i) = ((-1)**abs(m))*(Fact(l-m)/Fact(l+m))*Ps(el+1,el+1+em) end if @@ -111,4 +112,4 @@ real(db) function Y_lm(l,m,x,y,z)result(Ylm) end function Y_lm -End Module Spherical_Harmonics \ No newline at end of file +End Module Spherical_Harmonics From 71484087df229a65f74e0f1bf7bee0f97018aba7 Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Thu, 18 Sep 2025 16:03:35 +0100 Subject: [PATCH 09/17] Fix Fortitude rule C001 --- Code/ylm.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Code/ylm.f90 b/Code/ylm.f90 index c272454..5f6e948 100644 --- a/Code/ylm.f90 +++ b/Code/ylm.f90 @@ -8,7 +8,8 @@ !>@details !! It contains the definition \f$Y_{LM}\f$ in terms of Associate Legendre Polinomial (\f$P_{LM}\f$) Module Spherical_Harmonics - USE Params, ONLY: db,pi + USE Params, ONLY: db,pi + implicit none CONTAINS ! DESCRIPTION: Fact(n) !> @brief From 1b70c5b745a65c341fc25cb869da12230a48f01a Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Thu, 18 Sep 2025 16:05:54 +0100 Subject: [PATCH 10/17] Fix Fortitude rule C141 --- Code/twobody.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Code/twobody.f90 b/Code/twobody.f90 index c6abcbe..48f29a5 100644 --- a/Code/twobody.f90 +++ b/Code/twobody.f90 @@ -601,10 +601,10 @@ SUBROUTINE cut_system(full_analysis) center(:,ifrag)=center(:,ifrag)/mass(ifrag) END DO ! stop if this is not a separated case - IF(.NOT.istwobody) EXIT + IF(.NOT.istwobody) EXIT Iteration ! end iterations if the c.m. remains constant IF(MAX(MAXVAL(ABS(cent(1,:)-center(1,:))), & - MAXVAL(ABS(cent(3,:)-center(3,:))))<1.0d-05) EXIT + MAXVAL(ABS(cent(3,:)-center(3,:))))<1.0d-05) EXIT Iteration ENDDO Iteration ! record separation distance IF(istwobody) THEN From 68379e65caa9b9a6ffd8f30aa95b016a7e75cbed Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Thu, 18 Sep 2025 20:44:33 +0100 Subject: [PATCH 11/17] Fix Fortitude rule C121 --- Code/coulomb.f90 | 2 +- Code/dynamic.f90 | 22 ++++++++++++++-------- Code/energies.f90 | 9 ++++++--- Code/external.f90 | 2 +- Code/fourier.f90 | 2 +- Code/fragments.f90 | 10 ++++++---- Code/inout.f90 | 10 +++++++--- Code/levels.f90 | 5 +++-- Code/main3d.f90 | 12 ++++++------ Code/meanfield.f90 | 6 +++--- Code/moment.f90 | 6 ++++-- Code/pairs.f90 | 2 +- Code/parallel.f90 | 1 - Code/static.f90 | 18 ++++++++++++------ Code/twobody.f90 | 8 ++++---- Code/user.f90 | 3 +++ Code/user_sample.f90 | 5 ++--- 17 files changed, 74 insertions(+), 49 deletions(-) diff --git a/Code/coulomb.f90 b/Code/coulomb.f90 index a06bbc7..8acb562 100644 --- a/Code/coulomb.f90 +++ b/Code/coulomb.f90 @@ -111,7 +111,7 @@ 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 diff --git a/Code/dynamic.f90 b/Code/dynamic.f90 index 73c6e71..f1963ff 100644 --- a/Code/dynamic.f90 +++ b/Code/dynamic.f90 @@ -9,19 +9,25 @@ !!are equally used in the static calculation. !------------------------------------------------------------------------------ MODULE DYNAMIC - USE Params - USE Grids, ONLY: nx,ny,nz,wxyz - USE Densities - USE Levels - USE Energies - USE Moment + USE Params, ONLY: energiesfile, iter, momentafile, monopolesfile, mrest, & + octupolesfile, quadrupolesfile, spinfile, trestart, energiesfile, & + mplot, mprint, nof, printnow, diatriacontapolesfile, dipolesfile, & + hexadecapolesfile,hbc,db,extfieldfile,wflag,scratch,time + USE Grids, ONLY: nx,ny,nz,wxyz,x,y,z + USE Densities, ONLY: rho, tau, sdens, sodens, current + USE Levels, ONLY: wocc, npmin, psi, sp_norm, sp_energy, sp_orbital, sp_spin, & + sp_kinetic, nstmax,nstloc,isospin + USE Forces, ONLY:f + USE Energies, ONLY: total_angmom, tke, ehfls, ehf, ehf0, ehf1, ehf2, ehf3, & + ehfc, ehfint, orbital, spin + USE Moment, ONLY: pnr, pcm USE Twobody, ONLY: twobody_analysis,istwobody,roft,roft_old - USE Parallel + USE Parallel, ONLY: tmpi, mpi_nprocs,globalindex USE Meanfield, ONLY: skyrme, hpsi, spot USE Trivial, ONLY: overlap USE Inout, ONLY: write_wavefunctions,write_densities, plot_density, & sp_properties,start_protocol - USE External + USE External, ONLY: l_val, M_val IMPLICIT NONE SAVE INTEGER :: nt !< the number of the final time step to be diff --git a/Code/energies.f90 b/Code/energies.f90 index 34e08be..6d0e0de 100644 --- a/Code/energies.f90 +++ b/Code/energies.f90 @@ -69,9 +69,12 @@ !------------------------------------------------------------------------------ MODULE Energies USE Params, ONLY: db,tcoul - USE Forces - USE Densities - USE Levels + USE Forces, ONLY: f, b0, b0p, b1, b1p, b2, b2p, b3, b3p, b4, b4, b4p, slate, & + ipair + USE Densities, ONLY: rho, tau, sdens, current, sodens + USE Levels, ONLY: sp_kinetic, sp_energy, sp_efluct1, sp_efluct2, wocc, sp_orbital, & + sp_spin + USE Grids, ONLY: nx,ny,nz USE Pairs, ONLY: epair IMPLICIT NONE SAVE diff --git a/Code/external.f90 b/Code/external.f90 index ca254e2..cb829c4 100644 --- a/Code/external.f90 +++ b/Code/external.f90 @@ -42,7 +42,7 @@ MODULE External USE Grids, ONLY: nx,ny,nz,x,y,z,dx,dy,dz,wxyz USE Levels, ONLY: nstloc,isospin,charge_number,mass_number USE MEANFIELD, ONLY: upot - Use Spherical_Harmonics + Use Spherical_Harmonics, ONLY: y_lm IMPLICIT NONE INTEGER :: isoext=0 !< isospin behavior of the external field: \c isoext !! denotes the same action on protons and neutrons, diff --git a/Code/fourier.f90 b/Code/fourier.f90 index 6f57340..203f618 100644 --- a/Code/fourier.f90 +++ b/Code/fourier.f90 @@ -19,7 +19,7 @@ MODULE Fourier USE params, ONLY: db,wflag USE Grids, ONLY: nx,ny,nz - USE ISO_C_BINDING + USE ISO_C_BINDING, ONLY: c_long IMPLICIT NONE !>@name plans for full three-dimensional forward and backward transforms for both spin components. !>@{ diff --git a/Code/fragments.f90 b/Code/fragments.f90 index b8b818a..1d625c9 100644 --- a/Code/fragments.f90 +++ b/Code/fragments.f90 @@ -28,10 +28,12 @@ !------------------------------------------------------------------------------ MODULE Fragments - USE Params - USE Grids + USE Params, ONLY: db,mnof,wflag,e2,hbc,pi,scratch,wffile,nof,iter,tdynamic,time, & + trestart,tstatic + USE Grids, ONLY: x,y,z,nx,ny,nz,dx,dy,dz USE Forces, ONLY: f,nucleon_mass - USE Levels + USE Levels, ONLY: npmin, npsi, sp_efluct1, sp_kinetic, sp_norm, sp_parity, sp_energy, & + wocc, isospin, psi, nstmax, nprot, nneut, mass_number, charge_number IMPLICIT NONE SAVE PRIVATE @@ -280,7 +282,7 @@ END SUBROUTINE read_fragments !--------------------------------------------------------------------------- SUBROUTINE read_one_fragment(iff) USE Parallel, ONLY: node,mpi_myproc,localindex - USE Fourier + USE Fourier, ONLY: pbackward,pforward INTEGER,INTENT(IN) :: iff LOGICAL :: multifile INTEGER :: ipn diff --git a/Code/inout.f90 b/Code/inout.f90 index 6679fb7..d72d776 100644 --- a/Code/inout.f90 +++ b/Code/inout.f90 @@ -30,15 +30,19 @@ !!match the purpose of this module but are placed here for convenience. !------------------------------------------------------------------------------ MODULE Inout - USE Params + USE Params, ONLY: scratch, trestart, tfft, write_isospin, iter, nselect, & + tdynamic, time, scratch2, wffile, writeselect,db,wflag USE Parallel, ONLY: node,localindex,mpi_myproc - USE Grids + USE Grids, ONLY: x,y,z,der1x,der2x,der1y,der2y,der1z,der2z,wxyz,& + nx,ny,nz,dx,dy,dz USE Forces, ONLY:f USE Moment, ONLY: cm,cmtot USE Densities, ONLY: rho,tau,current,sdens,sodens USE Meanfield, ONLY: upot USE Coulomb, ONLY: wcoul - USE Levels + USE Levels, ONLY: sp_kinetic,sp_orbital,sp_parity,isospin,wocc,sp_norm, & + sp_energy,sp_efluct1,nstmax,npsi,nprot,nneut,npmin,charge_number,& + sp_spin,mass_number,psi IMPLICIT NONE PRIVATE :: write_one_density,write_vec_density CONTAINS diff --git a/Code/levels.f90 b/Code/levels.f90 index 7b17d42..4cb0b89 100644 --- a/Code/levels.f90 +++ b/Code/levels.f90 @@ -7,9 +7,10 @@ !!the pertinent arrays, allocating their storage and simple operations on them. !------------------------------------------------------------------------------ MODULE Levels - USE Params, ONLY: db,pi + USE Params, ONLY: db,pi, wflag USE Grids, ONLY: nx,ny,nz,dx,dy,dz - USE Fourier + USE Fourier, ONLY: pbackward, pforward, zbackward, zforward, ybackward, & + yforward, xbackward, xforward IMPLICIT NONE SAVE INTEGER :: nstmax ! Date: Thu, 18 Sep 2025 20:47:48 +0100 Subject: [PATCH 12/17] Fix Fortitude rule C081 --- Code/forces.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Code/forces.f90 b/Code/forces.f90 index 0b6d902..05a2eb8 100644 --- a/Code/forces.f90 +++ b/Code/forces.f90 @@ -130,11 +130,12 @@ SUBROUTINE read_force REAL(db) :: v0prot,v0neut,rho0pr,mixture INTEGER :: i LOGICAL :: predefined - LOGICAL :: turnoff_zpe=.FALSE. + LOGICAL :: turnoff_zpe ! read force definition NAMELIST /force/ name,pairing, & ex,zpe,h2m,t0,t1,t2,t3,t4,x0,x1,x2,x3,b4p,power, & ipair,v0prot,v0neut,rho0pr,mixture,turnoff_zpe + turnoff_zpe=.FALSE. ! mark force & pairing parameters as undefined h2m=-1.0; v0prot=-1.0; v0neut=-1.0; READ(5,force) From b502a2cfaddff7d730c80411a6f18d2de1a2e321 Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Thu, 18 Sep 2025 20:50:36 +0100 Subject: [PATCH 13/17] Fix Fortitude rule C072 --- Code/grids.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Code/grids.f90 b/Code/grids.f90 index 07f4ae7..6e84f5f 100644 --- a/Code/grids.f90 +++ b/Code/grids.f90 @@ -136,7 +136,7 @@ END SUBROUTINE init_grid !> REAL(db), array, returns matrix for damping. !--------------------------------------------------------------------------- SUBROUTINE init_coord(name,nv,dv,v,der1v,der2v,cdmpv) - CHARACTER(*) :: name + CHARACTER(len=1) :: name INTEGER :: nv REAL(db) :: dv INTENT(IN) :: name,nv,dv @@ -146,7 +146,7 @@ SUBROUTINE init_coord(name,nv,dv,v,der1v,der2v,cdmpv) ALLOCATE(v(nv),der1v(nv,nv),der2v(nv,nv),cdmpv(nv,nv)) v=[ ((i-1)*dv-0.5D0*FLOAT(nv-1)*dv,i=1,nv) ] IF(wflag) THEN - WRITE(*,'(1X,A,I3,A,F8.4,2(A,F8.4))') name // ' direction: ',nv, & + WRITE(*,'(1X,A,I3,A,F8.4,2(A,F8.4))') trim(name) // ' direction: ',nv, & ' points, spacing:',dv,' ranging from ',v(1),' to ',v(nv) ENDIF CALL sder(der1v,nv,dv) From cc205e082dac394c595adec16f93ffb367d109e2 Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Thu, 18 Sep 2025 20:54:38 +0100 Subject: [PATCH 14/17] Fix Fortitude rule PORT011, C071 and C002 --- Code/sequential.f90 | 2 +- Code/static.f90 | 9 ++++++--- Utils/Convert/Tdhf2Silo.f90 | 2 +- Utils/Overlap/overlap.f90 | 20 ++++++++++---------- 4 files changed, 18 insertions(+), 15 deletions(-) diff --git a/Code/sequential.f90 b/Code/sequential.f90 index cd6d58b..0918823 100644 --- a/Code/sequential.f90 +++ b/Code/sequential.f90 @@ -79,7 +79,7 @@ SUBROUTINE mpi_allreduce(rho,tmp_rho,length, & i_double_precision,sum, & comm_world,ierror) INTEGER :: ierror, comm_world, i_double_precision, length, sum - REAL(db), DIMENSION(*), INTENT(IN) :: rho,tmp_rho + REAL(db), DIMENSION(:), INTENT(IN) :: rho,tmp_rho STOP ' parallel calls inhibited ' RETURN END SUBROUTINE mpi_allreduce diff --git a/Code/static.f90 b/Code/static.f90 index 8b35eba..e7e9faa 100644 --- a/Code/static.f90 +++ b/Code/static.f90 @@ -526,10 +526,13 @@ SUBROUTINE diagstep(iq,nlin) SUBROUTINE zheevd( jobz, uplo, n, a, lda, w, work, & lwork, rwork, lrwork, iwork, liwork, info ) USE Params, ONLY: db + implicit none + INTEGER :: lda CHARACTER(1) :: jobz, uplo - INTEGER :: info, ldab, liwork, lrwork, lwork, n, iwork(*) - DOUBLE PRECISION :: rwork( * ), w( * ) - COMPLEX(8) :: a( lda, * ), work( * ) + INTEGER :: info, ldab, liwork, lrwork, lwork, n, iwork(:) + DOUBLE PRECISION :: rwork(:), w(:) + COMPLEX(db) :: a(:,:) + COMPLEX(db) :: work(:) INTENT(IN) :: jobz,uplo,n,lda,lwork,lrwork,liwork INTENT(INOUT) :: a INTENT(OUT) :: w,work,rwork,iwork,info diff --git a/Utils/Convert/Tdhf2Silo.f90 b/Utils/Convert/Tdhf2Silo.f90 index 6409934..b949914 100755 --- a/Utils/Convert/Tdhf2Silo.f90 +++ b/Utils/Convert/Tdhf2Silo.f90 @@ -108,7 +108,7 @@ PROGRAM Tdhf2Silo CONTAINS SUBROUTINE write_density(rname,data) CHARACTER(len=*),INTENT(IN) :: rname - REAL(8),INTENT(IN) :: data(*) + REAL(db),INTENT(IN) :: data(:) WRITE(*,*) 'Writing ',rname,' ',LEN(rname) ret=dbputqv1(dbid,rname,LEN(rname),'Mesh',4,data,dims,3, & DB_F77NULL,0,DB_DOUBLE,DB_NODECENT,op,iret) diff --git a/Utils/Overlap/overlap.f90 b/Utils/Overlap/overlap.f90 index e02c765..d908022 100755 --- a/Utils/Overlap/overlap.f90 +++ b/Utils/Overlap/overlap.f90 @@ -3,13 +3,13 @@ PROGRAM Analyze INTEGER :: nstmax,nneut,nprot,nx,ny,nz INTEGER :: iter,number(2),npsi(2) INTEGER :: nstmax2,nneut2,nprot2,nx2,ny2,nz2 - REAL(8) :: time,charge_number,mass_number,cm(3),cm2(3),dx,dy,dz, & + REAL(db) :: time,charge_number,mass_number,cm(3),cm2(3),dx,dy,dz, & wxyz,wxyz2 - CHARACTER(8) :: forcename,forcename2 - COMPLEX(8),ALLOCATABLE,DIMENSION(:,:,:,:,:) :: psi1,psi2 - COMPLEX(8) :: norm,on,op - COMPLEX(8),ALLOCATABLE :: ovn(:,:),ovp(:,:) - REAL(8) :: fn,fp,dnmax,dnmin,dpmax,dpmin,onmax,opmax + CHARACTER(db) :: forcename,forcename2 + COMPLEX(db),ALLOCATABLE,DIMENSION(:,:,:,:,:) :: psi1,psi2 + COMPLEX(db) :: norm,on,op + COMPLEX(db),ALLOCATABLE :: ovn(:,:),ovp(:,:) + REAL(db) :: fn,fp,dnmax,dnmin,dpmax,dpmin,onmax,opmax INTEGER :: i,j,m CHARACTER(40) :: file1,file2 DO @@ -106,14 +106,14 @@ PROGRAM Analyze END DO 1 STOP CONTAINS - COMPLEX(8) FUNCTION scalprod(a,b) + COMPLEX(db) FUNCTION scalprod(a,b) IMPLICIT NONE - COMPLEX(8),INTENT(IN),DIMENSION(:,:,:,:) :: a,b + COMPLEX(db),INTENT(IN),DIMENSION(:,:,:,:) :: a,b scalprod=SUM(CONJG(a)*b)*wxyz END FUNCTION scalprod - COMPLEX(8) function determinant(a,n) + COMPLEX(db) function determinant(a,n) INTEGER,INTENT(IN) :: n - COMPLEX(8) :: a(n,n),det(2),work(n) + COMPLEX(db) :: a(n,n),det(2),work(n) integer :: ipvt(n),info call zgefa(a,n,n,ipvt,info) if(info/=0) then From 9f4228944fc929bd12b4ccb3350e6b513b09493b Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Fri, 19 Sep 2025 13:17:14 +0100 Subject: [PATCH 15/17] Fix Fortitude rule C061 --- Code/fragments.f90 | 6 ++--- Code/grids.f90 | 23 ++++++++--------- Code/meanfield.f90 | 10 +++----- Code/moment.f90 | 3 ++- Code/sequential.f90 | 12 ++++----- Code/static.f90 | 24 +++++++++--------- Code/trivial.f90 | 52 +++++++++++++++++---------------------- Code/twobody.f90 | 6 ++--- Code/ylm.f90 | 15 ++++++----- Utils/Cuts/Cuts.f90 | 7 +++--- Utils/Overlap/overlap.f90 | 3 ++- 11 files changed, 76 insertions(+), 85 deletions(-) diff --git a/Code/fragments.f90 b/Code/fragments.f90 index 1d625c9..8741eae 100644 --- a/Code/fragments.f90 +++ b/Code/fragments.f90 @@ -433,10 +433,10 @@ END SUBROUTINE locate !> REAL(db), takes the shift !--------------------------------------------------------------------------- PURE SUBROUTINE phases(n,a,c) - INTEGER :: n,i,k + INTEGER :: i,k + INTEGER, intent(in) :: n COMPLEX(db),INTENT(OUT) :: a(n) - REAL(db) :: c - INTENT(IN) :: n,c + REAL(db), intent(in) :: c DO i=1,n k=i-1 IF(i>(n+1)/2) k=k-n diff --git a/Code/grids.f90 b/Code/grids.f90 index 6e84f5f..3695083 100644 --- a/Code/grids.f90 +++ b/Code/grids.f90 @@ -136,11 +136,10 @@ END SUBROUTINE init_grid !> REAL(db), array, returns matrix for damping. !--------------------------------------------------------------------------- SUBROUTINE init_coord(name,nv,dv,v,der1v,der2v,cdmpv) - CHARACTER(len=1) :: name - INTEGER :: nv - REAL(db) :: dv - INTENT(IN) :: name,nv,dv - REAL(db),POINTER :: v(:),der1v(:,:), & + CHARACTER(len=1), intent(in) :: name + INTEGER, intent(in) :: nv + REAL(db), intent(in) :: dv + REAL(db),POINTER, intent(out) :: v(:),der1v(:,:), & der2v(:,:),cdmpv(:,:) INTEGER :: i ALLOCATE(v(nv),der1v(nv,nv),der2v(nv,nv),cdmpv(nv,nv)) @@ -164,10 +163,9 @@ END SUBROUTINE init_coord !> REAL(db), takes the grid spacing. !--------------------------------------------------------------------------- PURE SUBROUTINE sder(der,nmax,d) - INTEGER :: nmax - REAL(db) :: d,der(:,:) - INTENT(IN) :: nmax,d - INTENT(OUT) :: der + INTEGER, intent(in) :: nmax + REAL(db), intent(in) :: d + REAL(db), intent(out) :: der(:,:) REAL(db) :: afac,sum INTEGER :: i,k,j,icn icn=(nmax+1)/2 @@ -195,10 +193,9 @@ END SUBROUTINE sder !> REAL(db), takes the grid spacing. !--------------------------------------------------------------------------- PURE SUBROUTINE sder2(der,nmax,d) - INTEGER :: nmax - REAL(db) :: d,der(1:nmax,1:nmax) - INTENT(IN) :: nmax,d - INTENT(OUT) :: der + INTEGER, intent(in) :: nmax + REAL(db), intent(in) :: d + REAL(db), intent(out) :: der(1:nmax,1:nmax) REAL(db) :: afac,sum INTEGER :: i,k,j,icn icn=(nmax+1)/2 diff --git a/Code/meanfield.f90 b/Code/meanfield.f90 index 3611fac..ea3a92e 100644 --- a/Code/meanfield.f90 +++ b/Code/meanfield.f90 @@ -318,12 +318,10 @@ END SUBROUTINE skyrme SUBROUTINE hpsi(iq,eshift,pinn,pout) USE Trivial, ONLY: cmulx, cmuly, cmulz USE Levels, ONLY: cdervx,cdervy,cdervz - INTEGER :: iq - REAL(db) :: eshift - COMPLEX(db),DIMENSION(:,:,:,:) :: pinn,pout - INTENT(IN) :: iq,eshift - INTENT(INOUT) :: pinn - INTENT(OUT) :: pout + INTEGER, intent(in) :: iq + REAL(db), intent(in) :: eshift + COMPLEX(db),DIMENSION(:,:,:,:) , intent(out) :: pinn + COMPLEX(db),DIMENSION(:,:,:,:) , intent(out) :: pout INTEGER :: is,ic REAL(db) :: sigis COMPLEX(db),ALLOCATABLE,DIMENSION(:,:,:,:) :: pswk,pswk2 diff --git a/Code/moment.f90 b/Code/moment.f90 index cc3cba1..56f2ab6 100644 --- a/Code/moment.f90 +++ b/Code/moment.f90 @@ -125,7 +125,8 @@ MODULE Moment !--------------------------------------------------------------------------- SUBROUTINE moments(L_val,M_val) USE Densities, ONLY: rho,current - INTEGER :: ix,iy,iz,iq,L_val,M_val,iii + INTEGER :: ix,iy,iz,iq,iii + INTEGER, intent(in) :: L_val,M_val REAL(db) :: xx(3),x2(3),vol,radius,eta REAL(db) :: qmat(3,3,2),qmtot(3,3) REAL(db) :: Mono(2),Quad(2),Oct(2),HexaDeca(2),DiaTriaConta(2),tmp,facn,facp diff --git a/Code/sequential.f90 b/Code/sequential.f90 index 0918823..f7bdb35 100644 --- a/Code/sequential.f90 +++ b/Code/sequential.f90 @@ -31,35 +31,35 @@ END SUBROUTINE init_all_mpi !************************************************************************ !> dummy function for the MPI routine SUBROUTINE mpi_init(ierror) - INTEGER :: ierror + INTEGER, intent(in) :: ierror STOP ' MPI_INIT: parallel calls inhibited ' RETURN END SUBROUTINE mpi_init !************************************************************************ !> dummy function for the MPI routine SUBROUTINE mpi_comm_size(comm_world,nprocs,ierror) - INTEGER :: ierror, nprocs, comm_world + INTEGER, intent(in) :: ierror, nprocs, comm_world STOP ' MPI_COMM_SIZE: parallel calls inhibited ' RETURN END SUBROUTINE mpi_comm_size !************************************************************************ !> dummy function for the MPI routine SUBROUTINE mpi_comm_rank(comm_world,myproc,ierror) - INTEGER :: ierror, myproc, comm_world + INTEGER, intent(in) :: ierror, myproc, comm_world STOP ' parallel calls inhibited ' RETURN END SUBROUTINE mpi_comm_rank !************************************************************************ !> dummy function for the MPI routine SUBROUTINE mpi_get_processor_name(processor_name,proc_namelen,ierror) - INTEGER :: ierror, processor_name, proc_namelen + INTEGER, intent(in) :: ierror, processor_name, proc_namelen STOP ' parallel calls inhibited ' RETURN END SUBROUTINE mpi_get_processor_name !************************************************************************ !> dummy function for the MPI routine SUBROUTINE mpi_barrier (comm_world, ierror) - INTEGER :: ierror, comm_world + INTEGER , intent(in) :: ierror, comm_world STOP ' parallel calls inhibited ' RETURN END SUBROUTINE mpi_barrier @@ -78,7 +78,7 @@ END SUBROUTINE associate_nodes SUBROUTINE mpi_allreduce(rho,tmp_rho,length, & i_double_precision,sum, & comm_world,ierror) - INTEGER :: ierror, comm_world, i_double_precision, length, sum + INTEGER, intent(in) :: ierror, comm_world, i_double_precision, length, sum REAL(db), DIMENSION(:), INTENT(IN) :: rho,tmp_rho STOP ' parallel calls inhibited ' RETURN diff --git a/Code/static.f90 b/Code/static.f90 index e7e9faa..1c4ed4d 100644 --- a/Code/static.f90 +++ b/Code/static.f90 @@ -427,10 +427,9 @@ SUBROUTINE grstep(nst,iq,spe,denerg,psin) ! !*********************************************************************** INTEGER,INTENT(IN) :: nst,iq - REAL(db) :: spe,denerg - COMPLEX(db) :: psin(:,:,:,:) - INTENT(OUT) :: denerg - INTENT(INOUT) :: spe,psin + REAL(db), intent(inout) :: spe + REAL(db), intent(out) :: denerg + COMPLEX(db), intent(inout) :: psin(:,:,:,:) REAL(db) :: x0act,esf,enrold,xnorm,xnormb,exph2,varh2 COMPLEX(db) :: ps1(nx,ny,nz,2),ps2(nx,ny,nz,2) INTEGER :: nst2 @@ -527,15 +526,14 @@ SUBROUTINE zheevd( jobz, uplo, n, a, lda, w, work, & lwork, rwork, lrwork, iwork, liwork, info ) USE Params, ONLY: db implicit none - INTEGER :: lda - CHARACTER(1) :: jobz, uplo - INTEGER :: info, ldab, liwork, lrwork, lwork, n, iwork(:) - DOUBLE PRECISION :: rwork(:), w(:) - COMPLEX(db) :: a(:,:) - COMPLEX(db) :: work(:) - INTENT(IN) :: jobz,uplo,n,lda,lwork,lrwork,liwork - INTENT(INOUT) :: a - INTENT(OUT) :: w,work,rwork,iwork,info + INTEGER, intent(in) :: lda + CHARACTER(1), intent(in) :: jobz, uplo + INTEGER :: ldab + INTEGER, intent(in) :: liwork, lrwork, lwork, n + integer, intent(out) :: info, iwork(:) + DOUBLE PRECISION, intent(out) :: rwork(:), w(:) + COMPLEX(db), intent(inout) :: a(:,:) + COMPLEX(db), intent(out) :: work(:) END SUBROUTINE zheevd END INTERFACE ! Step 1: copy matrix, then diagonalize diff --git a/Code/trivial.f90 b/Code/trivial.f90 index 97a9bed..a3979b4 100644 --- a/Code/trivial.f90 +++ b/Code/trivial.f90 @@ -33,11 +33,10 @@ MODULE Trivial !! input wave function. !--------------------------------------------------------------------------- PURE SUBROUTINE cmulx(xmat,pinn,pout,ifadd) - REAL(db) :: xmat(:,:) - COMPLEX(db) :: pinn(:,:,:,:),pout(:,:,:,:) - INTEGER :: ifadd - INTENT(IN) :: xmat,pinn,ifadd - INTENT(INOUT) :: pout + REAL(db), intent(in) :: xmat(:,:) + COMPLEX(db), intent(in) :: pinn(:,:,:,:) + COMPLEX(db), intent(inout) :: pout(:,:,:,:) + INTEGER, intent(in) :: ifadd INTEGER :: is,ix,iy,iz IF(ifadd==0) pout=0.0D0 FORALL(is=1:2,ix=1:nx,iy=1:ny,iz=1:nz) @@ -66,11 +65,10 @@ END SUBROUTINE cmulx !! input wave function. !--------------------------------------------------------------------------- PURE SUBROUTINE cmuly(ymat,pinn,pout,ifadd) - REAL(db) :: ymat(:,:) - COMPLEX(db) :: pinn(:,:,:,:),pout(:,:,:,:) - INTEGER :: ifadd - INTENT(IN) :: ymat,pinn,ifadd - INTENT(INOUT) :: pout + REAL(db), intent(in) :: ymat(:,:) + COMPLEX(db), intent(in) :: pinn(:,:,:,:) + COMPLEX(db), intent(inout) :: pout(:,:,:,:) + INTEGER, intent(in) :: ifadd INTEGER :: is,ix,iy,iz IF(ifadd==0) pout=0.0D0 FORALL(is=1:2,ix=1:nx,iy=1:ny,iz=1:nz) @@ -99,11 +97,10 @@ END SUBROUTINE cmuly !! input wave function. !--------------------------------------------------------------------------- PURE SUBROUTINE cmulz(zmat,pinn,pout,ifadd) - REAL(db) :: zmat(:,:) - COMPLEX(db) :: pinn(:,:,:,:),pout(:,:,:,:) - INTEGER :: ifadd - INTENT(IN) :: zmat,pinn,ifadd - INTENT(INOUT) :: pout + REAL(db), intent(in) :: zmat(:,:) + COMPLEX(db), intent(in) :: pinn(:,:,:,:) + COMPLEX(db), intent(inout) :: pout(:,:,:,:) + INTEGER, intent(in) :: ifadd INTEGER :: is,ix,iy,iz,izz IF(ifadd==0) pout=0.0D0 ! @@ -165,8 +162,8 @@ END FUNCTION rpsnorm !> COMPEX(db), array, takes \f$ psi_R \f$. !--------------------------------------------------------------------------- PURE FUNCTION overlap(pl,pr) RESULT(c) - COMPLEX(db) :: c,pl(:,:,:,:),pr(:,:,:,:) - INTENT(IN) :: pl,pr + COMPLEX(db), intent(in) :: pl(:,:,:,:),pr(:,:,:,:) + COMPLEX(db) :: c c=wxyz*SUM(CONJG(pl)*pr) END FUNCTION overlap !--------------------------------------------------------------------------- @@ -188,10 +185,9 @@ END FUNCTION overlap !! input wave function, if negative it is substracted. !--------------------------------------------------------------------------- PURE SUBROUTINE rmulx(xmat,finn,fout,ifadd) - INTEGER :: ifadd - REAL(db) :: xmat(:,:),finn(:,:,:),fout(:,:,:) - INTENT(IN) :: xmat,finn,ifadd - INTENT(INOUT) :: fout + INTEGER, intent(in) :: ifadd + REAL(db), intent(in) :: xmat(:,:),finn(:,:,:) + REAL(db), intent(inout) :: fout(:,:,:) INTEGER :: ix,iy,iz IF(ifadd==0) fout=0.D0 IF(ifadd>=0) THEN @@ -223,10 +219,9 @@ END SUBROUTINE rmulx !! input wave function, if negative it is substracted. !--------------------------------------------------------------------------- PURE SUBROUTINE rmuly(ymat,finn,fout,ifadd) - INTEGER :: ifadd - REAL(db) :: ymat(:,:),finn(:,:,:),fout(:,:,:) - INTENT(IN) :: ymat,finn,ifadd - INTENT(INOUT) :: fout + INTEGER, intent(in) :: ifadd + REAL(db), intent(in) :: ymat(:,:),finn(:,:,:) + REAL(db), intent(inout) :: fout(:,:,:) INTEGER :: ix,iy,iz IF(ifadd==0) fout=0.D0 IF(ifadd>=0) THEN @@ -258,10 +253,9 @@ END SUBROUTINE rmuly !! input wave function, if negative it is substracted. !--------------------------------------------------------------------------- PURE SUBROUTINE rmulz(zmat,finn,fout,ifadd) - INTEGER :: ifadd - REAL(db) :: zmat(:,:),finn(:,:,:),fout(:,:,:) - INTENT(IN) :: zmat,finn,ifadd - INTENT(INOUT) :: fout + INTEGER, intent(in) :: ifadd + REAL(db), intent(in) :: zmat(:,:),finn(:,:,:) + REAL(db), intent(inout) ::fout(:,:,:) INTEGER :: ix,iy,iz IF(ifadd==0) fout=0.D0 IF(ifadd>=0) THEN diff --git a/Code/twobody.f90 b/Code/twobody.f90 index 9b3b2c7..e6a6b38 100644 --- a/Code/twobody.f90 +++ b/Code/twobody.f90 @@ -234,9 +234,9 @@ END SUBROUTINE analyze !> REAL(db), solution vector. !--------------------------------------------------------------------------- SUBROUTINE gauss(a,b,c) - REAL(db) :: a(3,3),b(3),c(3),det - INTENT(IN) :: a,b - INTENT(OUT) :: c + REAL(db), intent(in) :: a(3,3),b(3) + REAL(db), intent(out) ::c(3) + REAL(db) :: det det=a(1,3)**2*a(2,2)-2.D0*a(1,2)*a(1,3)*a(2,3)+a(1,1)*a(2,3)**2+ & a(1,2)**2*a(3,3)-a(1,1)*a(2,2)*a(3,3) IF(ABS(det)<1.D-8) THEN diff --git a/Code/ylm.f90 b/Code/ylm.f90 index 5f6e948..a6c66f6 100644 --- a/Code/ylm.f90 +++ b/Code/ylm.f90 @@ -17,7 +17,8 @@ Module Spherical_Harmonics REAL(db) function Fact(n)result(Fct) IMPLICIT NONE - integer:: n,i + integer, intent(in):: n + integer :: i ! real(db)::Fct Fct=1.0d0 do i =1,n @@ -30,12 +31,13 @@ end function Fact REAL(db) function Fact2(n)result(Fct2) IMPLICIT NONE - integer :: n,i + integer, intent(in) :: n + integer :: i ! real(db)::Fct2 Fct2=1.0d0 - do while (n>0) - Fct2 = Fct2*n - n=n-2 + do i=0,n-1,2 + IF(i==n-1) exit + Fct2 = Fct2*(n-i) end do end function Fact2 @@ -46,7 +48,8 @@ real(db) function Plm(x,l,m)result(P_lm) IMPLICIT NONE integer, intent(in) :: l,m integer :: em_i,em_i_lower,em_i_lower2,el,em - real(db):: P00,P10,P11,x + real(db):: P00,P10,P11 + real(db), intent(in):: x real(db), ALLOCATABLE :: Ps(:,:) if (m>l)then diff --git a/Utils/Cuts/Cuts.f90 b/Utils/Cuts/Cuts.f90 index d350135..a20df16 100755 --- a/Utils/Cuts/Cuts.f90 +++ b/Utils/Cuts/Cuts.f90 @@ -40,10 +40,9 @@ PROGRAM Tdhf2Silo !******************************************************************* CONTAINS SUBROUTINE writeit(a,b,c,n,m,name) - REAL(db) :: a(:,:),b(:),c(:) - INTEGER :: n,m - CHARACTER(3) :: name - INTENT(IN) :: a,b,c,n,m,name + REAL(db), intent(in) :: a(:,:),b(:),c(:) + INTEGER, intent(in) :: n,m + CHARACTER(3), intent(in) :: name INTEGER :: i,j filename(7:9)=name WRITE(*,*) 'Writing ',filename diff --git a/Utils/Overlap/overlap.f90 b/Utils/Overlap/overlap.f90 index d908022..3531500 100755 --- a/Utils/Overlap/overlap.f90 +++ b/Utils/Overlap/overlap.f90 @@ -113,7 +113,8 @@ COMPLEX(db) FUNCTION scalprod(a,b) END FUNCTION scalprod COMPLEX(db) function determinant(a,n) INTEGER,INTENT(IN) :: n - COMPLEX(db) :: a(n,n),det(2),work(n) + COMPLEX(db), intent(out) :: a(n,n) + COMPLEX(db) :: det(2),work(n) integer :: ipvt(n),info call zgefa(a,n,n,ipvt,info) if(info/=0) then From 591096b89b243c247c509b04d27f606a894f748e Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Fri, 19 Sep 2025 14:19:26 +0100 Subject: [PATCH 16/17] Fix Fortitude rule C131 --- Code/coulomb.f90 | 2 +- Code/densities.f90 | 2 ++ Code/dynamic.f90 | 2 ++ Code/energies.f90 | 4 ++++ Code/external.f90 | 2 +- Code/forces.f90 | 3 +++ Code/fourier.f90 | 3 +++ Code/fragments.f90 | 1 - Code/grids.f90 | 4 +++- Code/inout.f90 | 4 +++- Code/levels.f90 | 4 ++++ Code/meanfield.f90 | 3 ++- Code/parallel.f90 | 3 +++ Code/params.f90 | 8 ++++++++ Code/sequential.f90 | 3 +++ Code/static.f90 | 2 ++ Code/trivial.f90 | 2 ++ Code/twobody.f90 | 2 ++ Code/user.f90 | 1 + Code/user_sample.f90 | 1 + Code/ylm.f90 | 2 ++ 21 files changed, 52 insertions(+), 6 deletions(-) diff --git a/Code/coulomb.f90 b/Code/coulomb.f90 index 8acb562..f4badbe 100644 --- a/Code/coulomb.f90 +++ b/Code/coulomb.f90 @@ -128,7 +128,7 @@ MODULE Coulomb COMPLEX(db),PRIVATE,ALLOCATABLE,SAVE :: q(:,:,:) !@name Scalar densities: !>These are dimensioned (nx,ny,nz,2), diff --git a/Code/dynamic.f90 b/Code/dynamic.f90 index f1963ff..3644b75 100644 --- a/Code/dynamic.f90 +++ b/Code/dynamic.f90 @@ -29,6 +29,8 @@ MODULE DYNAMIC sp_properties,start_protocol USE External, ONLY: l_val, M_val IMPLICIT NONE + PRIVATE + PUBLIC :: getin_dynamic,dynamichf SAVE INTEGER :: nt !< the number of the final time step to be !! calculated. In case of a restart this is smaller than the total number of time steps. diff --git a/Code/energies.f90 b/Code/energies.f90 index 6d0e0de..88d73b0 100644 --- a/Code/energies.f90 +++ b/Code/energies.f90 @@ -77,6 +77,10 @@ MODULE Energies USE Grids, ONLY: nx,ny,nz USE Pairs, ONLY: epair IMPLICIT NONE + PRIVATE + PUBLIC :: ipair, ehfprev, ehf, efluct2prev, efluct2, efluct1prev, & + efluct1,tke,total_angmom,spin,orbital,ehfls,ehfint,ehfc,ehf3,ehf2,ehf1, & + ehf0,ecorc,e3corr SAVE ! total energies calculated in subroutine "energy" and "hfenergy" REAL(db) :: ehft !< kinetic energy diff --git a/Code/external.f90 b/Code/external.f90 index cb829c4..1c2c2f4 100644 --- a/Code/external.f90 +++ b/Code/external.f90 @@ -80,7 +80,7 @@ MODULE External !! because there is no COM correction implemented for this, however, it works fine in the !! cases of big doubly magic nuclei like \f$^{208}Pb\f$. REAL(db) :: r_avg=0.0d0 !< The value of \f$ \sqrt{} \f$ used in the case of dipole boost. - + PRIVATE PUBLIC ::L_val,M_val,ampl_ext,isoext SAVE CONTAINS diff --git a/Code/forces.f90 b/Code/forces.f90 index 05a2eb8..3bfc1c4 100644 --- a/Code/forces.f90 +++ b/Code/forces.f90 @@ -21,6 +21,9 @@ MODULE Forces USE Params, ONLY: db,scratch,wflag, hbc, e2, pi IMPLICIT NONE + PRIVATE + PUBLIC :: h2ma,f,b0,b0p,b1,b1p,b2,b2p,b3,b3p,b4,b4p,slate,nucleon_mass, ipair,p, & + read_force SAVE !> This contains the parameters for pairing: TYPE Pairing diff --git a/Code/fourier.f90 b/Code/fourier.f90 index 203f618..72cb247 100644 --- a/Code/fourier.f90 +++ b/Code/fourier.f90 @@ -21,6 +21,9 @@ MODULE Fourier USE Grids, ONLY: nx,ny,nz USE ISO_C_BINDING, ONLY: c_long IMPLICIT NONE + PRIVATE + PUBLIC :: pbackward, pforward, zbackward, zforward, ybackward, yforward, & + xbackward, xforward !>@name plans for full three-dimensional forward and backward transforms for both spin components. !>@{ INTEGER(C_LONG),SAVE :: pforward,pbackward diff --git a/Code/fragments.f90 b/Code/fragments.f90 index 8741eae..13ace6d 100644 --- a/Code/fragments.f90 +++ b/Code/fragments.f90 @@ -83,7 +83,6 @@ MODULE Fragments INTEGER :: fnz(mnof) ! (nz,nz) and calculated in subroutine \c sder2. cdmpz ! (nz,nz) and calculated in subroutine \c setdmc. - PRIVATE :: init_coord, sder, sder2, setdmc, gauss + PRIVATE + PUBLIC :: nx,ny,nz,wxyz,dx,dy,dz,der1x,der1y,der1z,der2x,der2y,der2z, & + x,y,z,periodic,cdmpx,cdmpy,cdmpz,init_grid CONTAINS !--------------------------------------------------------------------------- ! DESCRIPTION: init_grid diff --git a/Code/inout.f90 b/Code/inout.f90 index d72d776..9ca84d7 100644 --- a/Code/inout.f90 +++ b/Code/inout.f90 @@ -44,7 +44,9 @@ MODULE Inout sp_energy,sp_efluct1,nstmax,npsi,nprot,nneut,npmin,charge_number,& sp_spin,mass_number,psi IMPLICIT NONE - PRIVATE :: write_one_density,write_vec_density + PRIVATE + PUBLIC :: write_wavefunctions,write_densities,plot_density,sp_properties, & + start_protocol CONTAINS !--------------------------------------------------------------------------- ! DESCRIPTION: write_wavefunctions diff --git a/Code/levels.f90 b/Code/levels.f90 index 4cb0b89..feea713 100644 --- a/Code/levels.f90 +++ b/Code/levels.f90 @@ -12,6 +12,10 @@ MODULE Levels USE Fourier, ONLY: pbackward, pforward, zbackward, zforward, ybackward, & yforward, xbackward, xforward IMPLICIT NONE + PRIVATE + PUBLIC :: cdervx,cdervy,cdervz,nstmax,npsi,nstloc,sp_kinetic,sp_orbital, & + sp_parity,isospin,wocc,sp_norm,sp_energy,sp_efluct1,nprot,nneut,npmin, & + charge_number,sp_spin,mass_number,psi,sp_efluct2,laplace,hmatr SAVE INTEGER :: nstmax ! @brief From 505f4a69a1153287dab6494e44cee09fa65008a7 Mon Sep 17 00:00:00 2001 From: mcbarton <150042563+mcbarton@users.noreply.github.com> Date: Fri, 19 Sep 2025 13:26:19 +0100 Subject: [PATCH 17/17] Ignore rule C003 Fortitude --- .github/workflows/fortitude.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/fortitude.yml b/.github/workflows/fortitude.yml index 15f0c3d..2a15c69 100644 --- a/.github/workflows/fortitude.yml +++ b/.github/workflows/fortitude.yml @@ -26,4 +26,4 @@ jobs: python-version: '3.13' - run: | python -m pip install fortitude-lint - fortitude check --output-format github + fortitude check --ignore C003 --output-format github