From 4f361a8dd2791d5a317a6e370d1157dd297fd333 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Wed, 8 Oct 2025 18:50:07 +0100 Subject: [PATCH] Update fortitude version and apply --preview fixes --- .github/workflows/format-lint.yml | 2 +- fortitude.toml | 2 ++ ruff.toml | 1 + src/compile_info.F90 | 2 +- src/force_bound.F90 | 4 ++-- src/force_h2o.F90 | 4 ++++ src/force_mm.F90 | 6 +++--- src/forces.F90 | 2 +- src/geom_analysis.F90 | 2 +- src/gle.F90 | 14 +++++++------- src/init.F90 | 2 ++ src/io.F90 | 4 ++-- src/landau_zener.F90 | 6 +++--- src/nosehoover.F90 | 32 +++++++++++++++---------------- src/plumed.F90 | 2 +- src/potentials.F90 | 10 +++++----- src/potentials_sh.F90 | 18 ++++++++--------- src/random.F90 | 28 +++++++++++++-------------- src/remd.F90 | 4 ++-- src/sh_integ.F90 | 6 +++--- src/surfacehop.F90 | 10 +++++----- src/transform.F90 | 12 ++++++------ src/utils.F90 | 8 ++++---- src/vinit.F90 | 8 ++++---- 24 files changed, 99 insertions(+), 90 deletions(-) create mode 100644 ruff.toml diff --git a/.github/workflows/format-lint.yml b/.github/workflows/format-lint.yml index 675876eb..0c6e39da 100644 --- a/.github/workflows/format-lint.yml +++ b/.github/workflows/format-lint.yml @@ -30,7 +30,7 @@ jobs: - uses: actions/checkout@v5 - name: Install fortitude linter - run: pip install fortitude-lint==0.7.1 + run: pip install fortitude-lint==0.7.5 - name: Lint check run: fortitude check src/ diff --git a/fortitude.toml b/fortitude.toml index a4e77794..0df73177 100644 --- a/fortitude.toml +++ b/fortitude.toml @@ -3,9 +3,11 @@ ignore = [ "line-too-long", "trailing-whitespace", # we should enable this, autofix available "incorrect-space-before-comment", # we should enable this, autofix available + "superfluous-semicolon", # same as above "use-all", # should use "only:" "assumed-size-character-intent", # TODO: We should fix all instances! "implicit-external-procedures", # requires Fortran 2018 + "bad-quote-string", # Use single or double quotes consistently ] exclude = [ "src/fftw3.F90", diff --git a/ruff.toml b/ruff.toml new file mode 100644 index 00000000..f11cf635 --- /dev/null +++ b/ruff.toml @@ -0,0 +1 @@ +line-length = 120 diff --git a/src/compile_info.F90 b/src/compile_info.F90 index 9423d5ab..b96d04db 100644 --- a/src/compile_info.F90 +++ b/src/compile_info.F90 @@ -4,7 +4,7 @@ ! and the current Git commit, which is passed in from Makefile. ! allow(procedure-not-in-module) ! fortitude linter subroutine print_compile_info() - use iso_fortran_env, only: compiler_version, compiler_options + use, intrinsic :: iso_fortran_env, only: compiler_version, compiler_options use mod_files, only: stdout character(len=*), parameter :: ABIN_VERSION = '1.2.0' diff --git a/src/force_bound.F90 b/src/force_bound.F90 index ab4aa047..e480fa20 100644 --- a/src/force_bound.F90 +++ b/src/force_bound.F90 @@ -41,7 +41,7 @@ subroutine sbc_init(x, y, z) ! calculation of the cluster radius, taken from the center of mass do iat = 1, natom r = (x(iat, iw) - xcm)**2 + (y(iat, iw) - ycm)**2 + (z(iat, iw) - zcm)**2 - r = dsqrt(r) + r = sqrt(r) if (r > rmax .and. names(iat) /= 'H') then rmax = r end if @@ -97,7 +97,7 @@ subroutine force_sbc(x, y, z, fx, fy, fz, walkmax) do iw = 1, walkmax do iat = 1, natom r = (x(iat, iw) - xcm)**2 + (y(iat, iw) - ycm)**2 + (z(iat, iw) - zcm)**2 - r = dsqrt(r) + r = sqrt(r) if (r > rb_sbc .and. names(iat) /= 'H') then frb = -kb_sbc * (r - rb_sbc) fx(iat, iw) = fx(iat, iw) + frb * (x(iat, iw) - xcm) / (r) diff --git a/src/force_h2o.F90 b/src/force_h2o.F90 index 733332e8..a22f79ac 100644 --- a/src/force_h2o.F90 +++ b/src/force_h2o.F90 @@ -192,6 +192,8 @@ subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads) y_new_forward(i, j) = y_new_forward(i, j) + DELTA case (3) z_new_forward(i, j) = z_new_forward(i, j) + DELTA + case default + call fatal_error(__FILE__, __LINE__, "Whoops, I should not be here!") end select ! Calculate the energy for the forward perturbed geometry @@ -209,6 +211,8 @@ subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads) fy(i, j) = -(Epot_delta(1) - Eclas_orig) / DELTA case (3) fz(i, j) = -(Epot_delta(1) - Eclas_orig) / DELTA + case default + call fatal_error(__FILE__, __LINE__, "Whoops, I should not be here!") end select end do end do diff --git a/src/force_mm.F90 b/src/force_mm.F90 index 3776d636..bc9821d4 100644 --- a/src/force_mm.F90 +++ b/src/force_mm.F90 @@ -181,7 +181,7 @@ subroutine initialize_LJ(rmin, eps, num_types, natom) if (LJcomb == 'LB') then ! WARNING: We expect rmin in angstroms! rij = 0.5D0 * (rmin(i) + rmin(j)) * ANG - epsij = dsqrt(eps(i) * eps(j)) + epsij = sqrt(eps(i) * eps(j)) end if Bij(i, j) = 2 * 6 * epsij * rij**6 Aij(i, j) = 12 * epsij * rij**12 @@ -211,7 +211,7 @@ subroutine force_mm(x, y, z, fx, fy, fz, eclas, walkmax) i = inames(iat1) j = inames(iat2) kLJ = ri3 * (ri3 * Aij(i, j) - Bij(i, j)) * ri - kC = q(i) * q(j) * dsqrt(ri3) + kC = q(i) * q(j) * sqrt(ri3) fx(iat1, iw) = fx(iat1, iw) + (kLJ + kC) * dx fx(iat2, iw) = fx(iat2, iw) - (kLJ + kC) * dx fy(iat1, iw) = fy(iat1, iw) + (kLJ + kC) * dy @@ -219,7 +219,7 @@ subroutine force_mm(x, y, z, fx, fy, fz, eclas, walkmax) fz(iat1, iw) = fz(iat1, iw) + (kLJ + kC) * dz fz(iat2, iw) = fz(iat2, iw) - (kLJ + kC) * dz eclas = eclas + ri3 * (ri3 * Aij(i, j) / 12 - Bij(i, j) / 6) - eclas = eclas + q(i) * q(j) / dsqrt(r) + eclas = eclas + q(i) * q(j) / sqrt(r) end do end do end do diff --git a/src/forces.F90 b/src/forces.F90 index 21006d0f..81a5eb1c 100644 --- a/src/forces.F90 +++ b/src/forces.F90 @@ -234,7 +234,7 @@ subroutine force_quantum(fx, fy, fz, x, y, z, amg, quantum_energy) ! Tuckerman normal modes Hamiltonian ! Not tested if (inormalmodes == 2) then - ak = NWALK * TEMP**2 * amg / dsqrt(nwalk * 1.0D0) + ak = NWALK * TEMP**2 * amg / sqrt(nwalk * 1.0D0) end if ! This is the energy coming from the Path Integral harmonic forces diff --git a/src/geom_analysis.F90 b/src/geom_analysis.F90 index 58feeac5..7408a06a 100644 --- a/src/geom_analysis.F90 +++ b/src/geom_analysis.F90 @@ -33,7 +33,7 @@ subroutine print_distances(x, y, z) ! If we have only one atom, take the distance from the origin, ! i.e. [0, 0, 0] if (natom == 1) then - r(idist) = dsqrt(x(1, iw)**2 + y(1, iw)**2 + z(1, iw)**2) + r(idist) = sqrt(x(1, iw)**2 + y(1, iw)**2 + z(1, iw)**2) else r(idist) = get_distance(x, y, z, dist1(idist), dist2(idist), iw) end if diff --git a/src/gle.F90 b/src/gle.F90 index e8a4bb83..3664e8ce 100644 --- a/src/gle.F90 +++ b/src/gle.F90 @@ -84,14 +84,14 @@ subroutine pile_init(dt, tau0) tau = tau0 / AUtoFS * 1000 gam = 1 / tau c1(1) = exp(-dt * gam) - c2(1) = dsqrt(1 - c1(1)**2) * dsqrt(temp * nwalk) + c2(1) = sqrt(1 - c1(1)**2) * sqrt(temp * nwalk) ! Now normal modes of bead necklace do iw = 2, nwalk omega = 2 * temp * nwalk * sin((iw - 1) * PI / nwalk) gam = 2 * omega c1(iw) = exp(-dt * gam) - c2(iw) = dsqrt(1 - c1(iw)**2) * dsqrt(temp * nwalk) + c2(iw) = sqrt(1 - c1(iw)**2) * sqrt(temp * nwalk) end do end subroutine pile_init @@ -116,9 +116,9 @@ subroutine pile_step(px, py, pz, m) ! TODO: implement global version according to equations 51-54 do iw = 1, nwalk do iat = 1, natom - px(iat, iw) = c1(iw) * px(iat, iw) + c2(iw) * ran(pom) * dsqrt(m(iat, iw)) - py(iat, iw) = c1(iw) * py(iat, iw) + c2(iw) * ran(pom + 1) * dsqrt(m(iat, iw)) - pz(iat, iw) = c1(iw) * pz(iat, iw) + c2(iw) * ran(pom + 2) * dsqrt(m(iat, iw)) + px(iat, iw) = c1(iw) * px(iat, iw) + c2(iw) * ran(pom) * sqrt(m(iat, iw)) + py(iat, iw) = c1(iw) * py(iat, iw) + c2(iw) * ran(pom + 1) * sqrt(m(iat, iw)) + pz(iat, iw) = c1(iw) * pz(iat, iw) + c2(iw) * ran(pom + 2) * sqrt(m(iat, iw)) pom = pom + 3 end do end do @@ -529,7 +529,7 @@ subroutine gle_propagate(p, T, S, mass, iw) do i = 1, ns + 1 call gautrg(ran, natom * 3) do j = 1, natom - sqm = dsqrt(mass(j, iw)) + sqm = sqrt(mass(j, iw)) p(j, i) = ran(j) * sqm p(j + natom, i) = ran(j + natom) * sqm p(j + natom * 2, i) = ran(j + natom * 2) * sqm @@ -610,7 +610,7 @@ subroutine cholesky(SST, S, n) end do do i = 1, n if (D(i, i) >= 0.0D0) then - D(i, i) = dsqrt(D(i, i)) + D(i, i) = sqrt(D(i, i)) else write (stderr, *) "WARNING: negative eigenvalue (", D(i, i), ")in LDL^T decomposition." D(i, i) = 0.0D0 diff --git a/src/init.F90 b/src/init.F90 index 19378f9e..d9c2bfe7 100644 --- a/src/init.F90 +++ b/src/init.F90 @@ -762,6 +762,8 @@ subroutine print_basic_info() write (stdout, *) ' Minimization ' case (5) write (stdout, *) ' Landau Zener MD ' + case default + call fatal_error(__FILE__, __LINE__, 'invalid ipimd in '//trim(chinput)) end select write (stdout, *) ' using potential: '//toupper(trim(pot)) diff --git a/src/io.F90 b/src/io.F90 index ffeda0d2..c3cf2c0c 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -45,7 +45,7 @@ subroutine print_dipoles(dipoles, iw, nstates) ! Most QM program print Dx, Dy, Dz components first. do ind = 1, 3 * nstates, 3 total_dip = dipoles(ind)**2 + dipoles(ind + 1)**2 + dipoles(ind + 2)**2 - total_dip = dsqrt(total_dip) + total_dip = sqrt(total_dip) write (UDIP, format2, advance="no") total_dip end do @@ -67,7 +67,7 @@ subroutine print_transdipoles(tdipoles, istate, ntdip) do ind = 1, 3 * ntdip, 3 total_tdip = tdipoles(ind)**2 + tdipoles(ind + 1)**2 + tdipoles(ind + 2)**2 - total_tdip = dsqrt(total_tdip) + total_tdip = sqrt(total_tdip) write (UTDIP, format2, advance="no") total_tdip end do diff --git a/src/landau_zener.F90 b/src/landau_zener.F90 index cd2dddab..83e6a983 100644 --- a/src/landau_zener.F90 +++ b/src/landau_zener.F90 @@ -255,7 +255,7 @@ subroutine lz_hop(x, y, z, vx, vy, vz, fxc, fyc, fzc, amt, dt, eclas, chpot) ! Three point minima of adiabatic splitting Zjk if ((en_diff(1) > en_diff(2)) .and. (en_diff(2) < en_diff(3)) .and. (it > 2)) then second_der = ((en_diff(3) - 2 * en_diff(2) + en_diff(1)) / dt**2) - prob(ist1) = exp(-PI / 2 * (dsqrt(en_diff(2)**3 / second_der))) + prob(ist1) = exp(-PI / 2 * (sqrt(en_diff(2)**3 / second_der))) write (fmt_in, '(I2.2)') ist write (fmt_out, '(I2.2)') ist1 write (stdout, *) "Three-point minimum (", trim(fmt_in), "->", trim(fmt_out), & @@ -364,7 +364,7 @@ subroutine lz_hop(x, y, z, vx, vy, vz, fxc, fyc, fzc, amt, dt, eclas, chpot) call force_clas(fxc, fyc, fzc, x, y, z, eclas, pot) !d) Simple velocity rescaling (https://doi.org/10.1063/1.4882073) - vel_rescale = dsqrt(1 - (dE / Ekin)) + vel_rescale = sqrt(1 - (dE / Ekin)) do iat = 1, natom vx(iat, 1) = vx(iat, 1) * vel_rescale vy(iat, 1) = vy(iat, 1) * vel_rescale @@ -398,7 +398,7 @@ subroutine lz_hop(x, y, z, vx, vy, vz, fxc, fyc, fzc, amt, dt, eclas, chpot) prob = 0.0D0 Ekin = ekin_v(vx, vy, vz) - molveloc = dsqrt(2.0D0 * Ekin / sum(am(1:natom))) + molveloc = sqrt(2.0D0 * Ekin / sum(am(1:natom))) call vranf(ran, 1) hop_rdnum = ran(1) icross = 0 diff --git a/src/nosehoover.F90 b/src/nosehoover.F90 index 336a461b..a86a6e60 100644 --- a/src/nosehoover.F90 +++ b/src/nosehoover.F90 @@ -301,9 +301,9 @@ subroutine initialize_nhc_momenta(temp) call gautrg(ran, natom * 3) ipom = 1 do iat = 1, natom - pnhx(iat, iw, inh) = ran(ipom) * dsqrt(temp * Qm(iw)) - pnhy(iat, iw, inh) = ran(ipom + 1) * dsqrt(temp * Qm(iw)) - pnhz(iat, iw, inh) = ran(ipom + 2) * dsqrt(temp * Qm(iw)) + pnhx(iat, iw, inh) = ran(ipom) * sqrt(temp * Qm(iw)) + pnhy(iat, iw, inh) = ran(ipom + 1) * sqrt(temp * Qm(iw)) + pnhz(iat, iw, inh) = ran(ipom + 2) * sqrt(temp * Qm(iw)) end do ipom = ipom + 3 end do @@ -316,7 +316,7 @@ subroutine initialize_nhc_momenta(temp) ! +1 if nmolt=1, gautrg needs array at least of length=2 call gautrg(ran, nmolt + 1) do imol = 1, nmolt - pnhx(imol, iw, inh) = ran(imol) * dsqrt(temp * ms(imol, inh)) + pnhx(imol, iw, inh) = ran(imol) * sqrt(temp * ms(imol, inh)) end do end do end do @@ -482,13 +482,13 @@ subroutine shiftNHC_yosh(px, py, pz, amt, dt) pnhx(imol, iw, nchain) = pnhx(imol, iw, nchain) + G(nchain) * wdt2 ! UPDATE THERMOSTAT VELOCITIES do inh = 1, nchain - 1 - AA = dexp(-wdt4 * pnhx(imol, iw, nchain - inh + 1) / ms(imol, nchain - inh + 1)) + AA = exp(-wdt4 * pnhx(imol, iw, nchain - inh + 1) / ms(imol, nchain - inh + 1)) pnhx(imol, iw, nchain - inh) = pnhx(imol, iw, nchain - inh) * AA * AA + & & wdt2 * G(nchain - inh) * AA end do ! UPDATE PARTICLE VELOCITIES - AA = dexp(-wdt * pnhx(imol, iw, 1) / ms(imol, 1)) + AA = exp(-wdt * pnhx(imol, iw, 1) / ms(imol, 1)) pscale = pscale * AA ! UPDATE FORCES G(1) = pscale * pscale * ekin2 - nf * temp @@ -501,7 +501,7 @@ subroutine shiftNHC_yosh(px, py, pz, amt, dt) ! UPDATE THERMOSTAT VELOCITIES do inh = 1, nchain - 1 - AA = dexp(-wdt4 * pnhx(imol, iw, inh + 1) / ms(imol, inh + 1)) + AA = exp(-wdt4 * pnhx(imol, iw, inh + 1) / ms(imol, inh + 1)) pnhx(imol, iw, inh) = pnhx(imol, iw, inh) * AA * AA + wdt2 * G(inh) * AA G(inh + 1) = pnhx(imol, iw, inh) * pnhx(imol, iw, inh) / ms(imol, inh) - temp end do @@ -564,23 +564,23 @@ subroutine shiftNHC_yosh_mass(px, py, pz, amt, dt) ! UPDATE THERMOSTAT VELOCITIES do inh = 1, nchain - 1 - AA = dexp(-wdt4 * pnhx(iat, iw, nchain - inh + 1) / Qm(iw)) + AA = exp(-wdt4 * pnhx(iat, iw, nchain - inh + 1) / Qm(iw)) pnhx(iat, iw, nchain - inh) = pnhx(iat, iw, nchain - inh) * AA * AA + & & wdt2 * Gx(nchain - inh) * AA - AA = dexp(-wdt4 * pnhy(iat, iw, nchain - inh + 1) / Qm(iw)) + AA = exp(-wdt4 * pnhy(iat, iw, nchain - inh + 1) / Qm(iw)) pnhy(iat, iw, nchain - inh) = pnhy(iat, iw, nchain - inh) * AA * AA + & & wdt2 * Gy(nchain - inh) * AA - AA = dexp(-wdt4 * pnhz(iat, iw, nchain - inh + 1) / Qm(iw)) + AA = exp(-wdt4 * pnhz(iat, iw, nchain - inh + 1) / Qm(iw)) pnhz(iat, iw, nchain - inh) = pnhz(iat, iw, nchain - inh) * AA * AA + & & wdt2 * Gz(nchain - inh) * AA end do ! UPDATE PARTICLE VELOCITIES - AA = dexp(-wdt * pnhx(iat, iw, 1) / Qm(iw)) + AA = exp(-wdt * pnhx(iat, iw, 1) / Qm(iw)) px(iat, iw) = px(iat, iw) * AA - AA = dexp(-wdt * pnhy(iat, iw, 1) / Qm(iw)) + AA = exp(-wdt * pnhy(iat, iw, 1) / Qm(iw)) py(iat, iw) = py(iat, iw) * AA - AA = dexp(-wdt * pnhz(iat, iw, 1) / Qm(iw)) + AA = exp(-wdt * pnhz(iat, iw, 1) / Qm(iw)) pz(iat, iw) = pz(iat, iw) * AA ! UPDATE FORCES @@ -597,13 +597,13 @@ subroutine shiftNHC_yosh_mass(px, py, pz, amt, dt) ! UPDATE THERMOSTAT VELOCITIES do inh = 1, nchain - 1 - AA = dexp(-wdt4 * pnhx(iat, iw, inh + 1) / Qm(iw)) + AA = exp(-wdt4 * pnhx(iat, iw, inh + 1) / Qm(iw)) pnhx(iat, iw, inh) = pnhx(iat, iw, inh) * AA * AA + wdt2 * Gx(inh) * AA Gx(inh + 1) = pnhx(iat, iw, inh)**2 / Qm(iw) - temp - AA = dexp(-wdt4 * pnhy(iat, iw, inh + 1) / Qm(iw)) + AA = exp(-wdt4 * pnhy(iat, iw, inh + 1) / Qm(iw)) pnhy(iat, iw, inh) = pnhy(iat, iw, inh) * AA * AA + wdt2 * Gy(inh) * AA Gy(inh + 1) = pnhy(iat, iw, inh)**2 / Qm(iw) - temp - AA = dexp(-wdt4 * pnhz(iat, iw, inh + 1) / Qm(iw)) + AA = exp(-wdt4 * pnhz(iat, iw, inh + 1) / Qm(iw)) pnhz(iat, iw, inh) = pnhz(iat, iw, inh) * AA * AA + wdt2 * Gz(inh) * AA Gz(inh + 1) = pnhz(iat, iw, inh)**2 / Qm(iw) - temp end do diff --git a/src/plumed.F90 b/src/plumed.F90 index 602a602b..362f830a 100644 --- a/src/plumed.F90 +++ b/src/plumed.F90 @@ -5,7 +5,7 @@ ! PLUMED is statically linked with ABIN, see Makefile. ! To install PLUMED libraries, see 'dev_scripts/install_plumed.sh' module mod_plumed - use iso_c_binding, only: C_NULL_PTR + use, intrinsic :: iso_c_binding, only: C_NULL_PTR use mod_const, only: DP #ifndef USE_PLUMED use mod_error, only: not_compiled_with diff --git a/src/potentials.F90 b/src/potentials.F90 index e4d018b7..c81ef65c 100644 --- a/src/potentials.F90 +++ b/src/potentials.F90 @@ -231,7 +231,7 @@ subroutine force_harmonic_rotor(x, y, z, fx, fy, fz, eclas, nwalk) dy = y(2, i) - y(1, i) dz = z(2, i) - z(1, i) r = dx**2 + dy**2 + dz**2 - r = dsqrt(r) + r = sqrt(r) fac = hrot%k * (r - hrot%r0) / r fx(1, i) = fac * dx fx(2, i) = -fx(1, i) @@ -258,7 +258,7 @@ subroutine hessian_harmonic_rotor(x, y, z, nwalk, hess) dy = y(2, i) - y(1, i) dz = z(2, i) - z(1, i) r = dx**2 + dy**2 + dz**2 - r = dsqrt(r) + r = sqrt(r) fac = hrot%k * (r - hrot%r0) / r hess(1, 1, i) = (hrot%k * dx**2 / r**2 - fac * dx**2 / r**2 + fac) / nwalk hess(2, 2, i) = (hrot%k * dy**2 / r**2 - fac * dy**2 / r**2 + fac) / nwalk @@ -313,7 +313,7 @@ subroutine morse_init(natom, k_morse, r0, d0) call real_positive(d0, 'd0_morse') call real_positive(r0, 'r0_morse') - a = dsqrt(k_morse / 2.0D0 / d0) + a = sqrt(k_morse / 2.0D0 / d0) morse = morse_params(d0=d0, a=a, r0=r0) end subroutine morse_init @@ -332,7 +332,7 @@ subroutine force_morse(x, y, z, fx, fy, fz, eclas, nwalk) dy = y(2, i) - y(1, i) dz = z(2, i) - z(1, i) r = dx**2 + dy**2 + dz**2 - r = dsqrt(r) + r = sqrt(r) ex = exp(-morse%a * (r - morse%r0)) fac = 2 * morse%a * ex * morse%d0 * (1 - ex) / r fx(1, i) = fac * dx @@ -364,7 +364,7 @@ subroutine hessian_morse(x, y, z, nwalk, hess) dy = y(2, i) - y(1, i) dz = z(2, i) - z(1, i) r = dx**2 + dy**2 + dz**2 - r = dsqrt(r) + r = sqrt(r) ex = exp(-a * (r - morse%r0)) fac = 2 * a * ex * morse%d0 * (1 - ex) / r fac2 = 2 * morse%d0 * a**2 * ex**2 / r**2 diff --git a/src/potentials_sh.F90 b/src/potentials_sh.F90 index 93fa1595..17f3dd7f 100644 --- a/src/potentials_sh.F90 +++ b/src/potentials_sh.F90 @@ -91,7 +91,7 @@ subroutine force_nai(x, y, z, fx, fy, fz, eclas) dy = y(2, 1) - y(1, 1) dz = z(2, 1) - z(1, 1) r = dx**2 + dy**2 + dz**2 - r = dsqrt(r) + r = sqrt(r) ! normalized distance dx = dx / r @@ -102,12 +102,12 @@ subroutine force_nai(x, y, z, fx, fy, fz, eclas) r = r / ANG ! calculating diabatic hamiltonian - VX = (nai%a2 + (nai%b2 / r)**8) * dexp(-r / nai%rho) - nai%e2 / r - nai%e2 * (nai%lp + nai%lm) / 2 / r**4 - & + VX = (nai%a2 + (nai%b2 / r)**8) * exp(-r / nai%rho) - nai%e2 / r - nai%e2 * (nai%lp + nai%lm) / 2 / r**4 - & nai%c2 / r**6 - 2 * nai%e2 * nai%lm * nai%lp / r**7 + nai%de0 - VA = nai%a1 * dexp(-nai%beta1 * (r - nai%r0)) - VXA = nai%a12 * dexp(-nai%beta12 * (r - nai%rx)**2) + VA = nai%a1 * exp(-nai%beta1 * (r - nai%r0)) + VXA = nai%a12 * exp(-nai%beta12 * (r - nai%rx)**2) ! calculating derivatives of diabatic hamiltonian - dVX = -(nai%a2 + (nai%b2 / r)**8) * dexp(-r / nai%rho) / nai%rho - 8.0D0 * nai%b2**8 / r**9 * dexp(-r / nai%rho) + & + dVX = -(nai%a2 + (nai%b2 / r)**8) * exp(-r / nai%rho) / nai%rho - 8.0D0 * nai%b2**8 / r**9 * exp(-r / nai%rho) + & 14.0D0 * nai%e2 * nai%lm * nai%lp / r**8 + 6.0D0 * nai%c2 / r**7 + 2.0D0 * nai%e2 * (nai%lp + nai%lm) / r**5 + & nai%e2 / r**2 dVA = -nai%beta1 * VA @@ -121,15 +121,15 @@ subroutine force_nai(x, y, z, fx, fy, fz, eclas) dVXA = dVXA / AUTOEV / ANG ! adiabatic potentials - E1 = (VA + VX) / 2.0D0 - dsqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0) / 2.0D0 - E2 = (VA + VX) / 2.0D0 + dsqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0) / 2.0D0 + E1 = (VA + VX) / 2.0D0 - sqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0) / 2.0D0 + E2 = (VA + VX) / 2.0D0 + sqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0) / 2.0D0 ! nonadiabatic coupling vector in the reduced system d12 = -(VXA * (dVA - dVX) + (-VA + VX) * dVXA) / (VA**2.0D0 - 2.0D0 * VA * VX + VX**2.0D0 + 4.0D0 * VXA**2.0D0) ! derivatives of energies dE1 = (dVA + dVX) / 2.0D0 - (2.0D0 * (VA - VX) * (dVA - dVX) + 8.0D0 * VXA * dVXA) / & - (4.0D0 * dsqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0)) + (4.0D0 * sqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0)) dE2 = (dVA + dVX) / 2.0D0 + (2.0D0 * (VA - VX) * (dVA - dVX) + 8.0D0 * VXA * dVXA) / & - (4.0D0 * dsqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0)) + (4.0D0 * sqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0)) ! saving electronic energies for SH en_array(1) = E1 diff --git a/src/random.F90 b/src/random.F90 index 8712ba3e..7ca7d57a 100644 --- a/src/random.F90 +++ b/src/random.F90 @@ -133,7 +133,7 @@ subroutine gautrg(gran, nran, iseed) -0.2017619095413939D-06, 0.2948964053761139D-05, & -0.1051448397925916D-06/ - if(isave.lt.0) then + if(isave<0) then isave=0 tiny = r1mach()**2 pi4=atan(one) @@ -149,9 +149,9 @@ subroutine gautrg(gran, nran, iseed) call vranf(gran, 0, iseed) end if - if(nran.gt.0) then + if(nran>0) then newran=nran-isave - if(isave.eq.1) gran(1)=gsave + if(isave==1) gran(1)=gsave call vranf(gran(isave+1), newran) do 100 i=1,newran-1,2 fac = sqrt(twom*log(gran(isave+i)+tiny)) @@ -174,7 +174,7 @@ subroutine gautrg(gran, nran, iseed) gran(isave+i+1)=fac*cxi 100 continue - if(mod(newran,2).eq.1) then + if(mod(newran,2)==1) then call vranf(xran, 1) fac=sqrt(twom*log(gran(nran)+tiny)) trig=twopi*xran(1) @@ -368,7 +368,7 @@ subroutine xuinit(y, np, nq, mode, nexec, iseed, init, last) ! set table to zero and exercise the bit generator a little - ix = iabs(iseed) + ix = abs(iseed) do i = 1, np x(i) = zero @@ -428,7 +428,7 @@ subroutine xuwarm(y, np, nq, mode, nexec) integer, intent(in) :: mode, nexec, np, nq integer :: i, k - if (nq.ge.np.or.np.eq.0.or.nq.eq.0) then + if (nq>=np.or.np==0.or.nq==0) then call fatal_error(__FILE__, __LINE__, & & 'Illegal table parameter(s) in xuwarm') return @@ -437,26 +437,26 @@ subroutine xuwarm(y, np, nq, mode, nexec) ! exercise the generator for nexec rounds of np prn's ! separate sections for subtractive or additive version - if(mode.le.0) then + if(mode<=0) then do k=1,nexec do i=1,nq y(i)=y(i)-y(i+np-nq) - if(y(i).lt.zero) y(i)=y(i)+one + if(y(i)=one) y(i)=y(i)-one end do do i=nq+1,np y(i)=y(i)+y(i-nq) - if(y(i).ge.one) y(i)=y(i)-one + if(y(i)>=one) y(i)=y(i)-one end do end do end if @@ -478,12 +478,12 @@ function R1MACH() ! THIS IS DEFINED AS THE SMALLEST POSITIVE MACHINE NUMBER ! U SUCH THAT 1.0 + U .NE. 1.0E0 ! --------------------------------------------------------------------- - IF(ICALL.EQ.0) THEN + IF(ICALL==0) THEN ICALL=1 U = ONE 10 U = U*HALF COMP = ONE + U - IF(COMP .NE. ONE) GO TO 10 + IF(COMP /= ONE) GOTO 10 EPS = U*TWO END IF R1MACH = EPS @@ -640,7 +640,7 @@ integer function get_random_seed() result(seed) write (stdout, *) 'Getting random seed from /dev/urandom' read (un) seed close (un) - seed = iabs(seed) + seed = abs(seed) else ! Fallback to XOR:ing the current time and pid. The PID is ! useful in case one launches multiple instances of the same diff --git a/src/remd.F90 b/src/remd.F90 index fa9d2db7..b6b8284f 100644 --- a/src/remd.F90 +++ b/src/remd.F90 @@ -175,9 +175,9 @@ subroutine swap_replicas(x, y, z, px, py, pz, fxc, fyc, fzc, rank1, rank2) ! Momenta are scaled according to the new temperature ! p_new = p_old * dsqrt(T_new/T_old) if (my_rank == rank1) then - scal = dsqrt(temp_list(my_rank + 1) / temp_list(my_rank + 2)) + scal = sqrt(temp_list(my_rank + 1) / temp_list(my_rank + 2)) else if (my_rank == rank2) then - scal = dsqrt(temp_list(my_rank + 1) / temp_list(my_rank)) + scal = sqrt(temp_list(my_rank + 1) / temp_list(my_rank)) end if if (my_rank == rank1) then diff --git a/src/sh_integ.F90 b/src/sh_integ.F90 index c932ae49..0bcacc87 100644 --- a/src/sh_integ.F90 +++ b/src/sh_integ.F90 @@ -371,14 +371,14 @@ subroutine sh_decoherence_correction(potential_energies, alpha, kinetic_energy, delta_e = abs(potential_energies(ist1) - potential_energies(current_state)) tau = (1.0D0 + alpha / kinetic_energy) / delta_e - scaling_factor = dexp(-dtp / tau) + scaling_factor = exp(-dtp / tau) ! WARNING: The Eq. 17 in the Persico paper is wrong, ! the scaling factor should be damping populations, NOT the WF coefficients. ! In previous ABIN versions (and in some other SH packages), this was incorrect. ! To reproduce older data, the user can set correct_decoherence=.false. in namelist sh. if (correct_decoherence) then - scaling_factor = dsqrt(scaling_factor) + scaling_factor = sqrt(scaling_factor) end if cel_re(ist1) = cel_re(ist1) * scaling_factor @@ -403,7 +403,7 @@ subroutine sh_decoherence_correction(potential_energies, alpha, kinetic_energy, call fatal_error(__FILE__, __LINE__, 'Invalid decoherence renormalization.') end if - renormalization_factor = dsqrt(renormalization_factor) + renormalization_factor = sqrt(renormalization_factor) cel_re(current_state) = cel_re(current_state) * renormalization_factor cel_im(current_state) = cel_im(current_state) * renormalization_factor diff --git a/src/surfacehop.F90 b/src/surfacehop.F90 index 921c9aeb..10b14e80 100644 --- a/src/surfacehop.F90 +++ b/src/surfacehop.F90 @@ -616,7 +616,7 @@ subroutine calc_baeckan(dt) de2dt2 = (2.0D0 * de(1) - 5.0D0 * de(2) + 4.0D0 * de(3) - de(4)) / dt**2 argument = de2dt2 / de(1) if (argument > 0.0D0) then - sigma_ba(ist2, ist1) = dsqrt(argument) / 2.0D0 + sigma_ba(ist2, ist1) = sqrt(argument) / 2.0D0 end if sigma_ba(ist1, ist2) = -sigma_ba(ist2, ist1) end do @@ -1044,11 +1044,11 @@ subroutine try_hop_nacme_rescale(vx, vy, vz, instate, outstate, eclas) ! Rescaling the velocities if (b_temp < 0) then - g_temp = (b_temp + dsqrt(c_temp)) / 2.0D0 / a_temp + g_temp = (b_temp + sqrt(c_temp)) / 2.0D0 / a_temp end if if (b_temp >= 0) then - g_temp = (b_temp - dsqrt(c_temp)) / 2.0D0 / a_temp + g_temp = (b_temp - sqrt(c_temp)) / 2.0D0 / a_temp end if iw = 1 @@ -1169,7 +1169,7 @@ subroutine try_hop_simple_rescale(vx, vy, vz, instate, outstate, eclas) if (ekin >= de) then - alfa = dsqrt(1 - de / ekin) + alfa = sqrt(1 - de / ekin) vx = alfa * vx vy = alfa * vy @@ -1296,7 +1296,7 @@ integer function check_CIVector(CIvecs, CIvecs_old, ci_len, nstates) do i = 1, ci_len cidotprod(ist1) = cidotprod(ist1) + CIvecs(i, ist1) * CIvecs_old(i, ist1) end do - if (cidotprod(ist1) < 1 / dsqrt(2.0D0)) then + if (cidotprod(ist1) < 1 / sqrt(2.0D0)) then write (*, *) "Warning: cidotprod too low." end if end do diff --git a/src/transform.F90 b/src/transform.F90 index 85c76a9e..edbb1be6 100644 --- a/src/transform.F90 +++ b/src/transform.F90 @@ -320,8 +320,8 @@ subroutine UtoX(x, y, z, transx, transy, transz) dnwalk = 1.0D0 fac = 1.0D0 if (inormalmodes == 1) then - dnwalk = dsqrt(1.0D0 * nwalk) - fac = dsqrt(2.0D0) + dnwalk = sqrt(1.0D0 * nwalk) + fac = sqrt(2.0D0) end if nmodes = nwalk / 2 @@ -386,8 +386,8 @@ subroutine XtoU(x, y, z, transx, transy, transz) dnwalk = 1.0D0 fac = 1.0D0 if (inormalmodes == 1) then - dnwalk = dsqrt(1.0D0 * nwalk) - fac = dsqrt(2.0D0) + dnwalk = sqrt(1.0D0 * nwalk) + fac = sqrt(2.0D0) end if nmodes = nwalk / 2 @@ -453,7 +453,7 @@ subroutine equant_cart(x, y, z, equant) omega_n = NWALK * TEMP if (inormalmodes == 2) then - omega_n = omega_n * dsqrt(NWALK * 1.D0) + omega_n = omega_n * sqrt(NWALK * 1.D0) end if equantx = 0.0D0 equanty = 0.0D0 @@ -497,7 +497,7 @@ subroutine equant_nm(x, y, z, equant) omega_n = NWALK * TEMP ! Tuckerman Hamiltonian if (inormalmodes == 2) then - omega_n = omega_n * dsqrt(NWALK * 1.D0) + omega_n = omega_n * sqrt(NWALK * 1.D0) end if equantx = 0.0D0 equanty = 0.0D0 diff --git a/src/utils.F90 b/src/utils.F90 index b3e46280..ce96ece0 100644 --- a/src/utils.F90 +++ b/src/utils.F90 @@ -22,7 +22,7 @@ real(DP) function get_distance(x, y, z, at1, at2, iw) result(r) r = (x(at1, iw) - x(at2, iw))**2 r = r + (y(at1, iw) - y(at2, iw))**2 r = r + (z(at1, iw) - z(at2, iw))**2 - r = dsqrt(r) + r = sqrt(r) end function get_distance real(DP) function get_angle(x, y, z, at1, at2, at3, iw) result(angle) @@ -46,7 +46,7 @@ real(DP) function get_angle(x, y, z, at1, at2, at3, iw) result(angle) vec2y = y(at3, iw) - y(at2, iw) vec2z = z(at3, iw) - z(at2, iw) angle = 180.0D0 / PI * acos((vec1x * vec2x + vec1y * vec2y + vec1z * vec2z) / & - & (dsqrt(vec1x**2 + vec1y**2 + vec1z**2) * dsqrt(vec2x**2 + vec2y**2 + vec2z**2))) + & (sqrt(vec1x**2 + vec1y**2 + vec1z**2) * sqrt(vec2x**2 + vec2y**2 + vec2z**2))) end function get_angle real(DP) function get_dihedral(x, y, z, at1, at2, at3, at4, iw, shiftdih) @@ -81,7 +81,7 @@ real(DP) function get_dihedral(x, y, z, at1, at2, at3, at4, iw, shiftdih) ! TODO: Add error handling for malformed dihedral angles to prevent division by zero get_dihedral = 180.0D0 / PI * acos( & & (norm1x * norm2x + norm1y * norm2y + norm1z * norm2z) / & - & (dsqrt(norm1x**2 + norm1y**2 + norm1z**2) * dsqrt(norm2x**2 + norm2y**2 + norm2z**2)) & + & (sqrt(norm1x**2 + norm1y**2 + norm1z**2) * sqrt(norm2x**2 + norm2y**2 + norm2z**2)) & & ) if (sign > 0) get_dihedral = shiftdih - get_dihedral @@ -109,7 +109,7 @@ end function sanitize_string ! Convert FORTRAN string to zero-terminated C string ! and remove any leading and trailing spaces. function c_string(string) - use iso_c_binding, only: C_CHAR, C_NULL_CHAR + use, intrinsic :: iso_c_binding, only: C_CHAR, C_NULL_CHAR character(kind=C_CHAR, len=*), intent(in) :: string character(kind=C_CHAR, len=len(string) + 1) :: c_string c_string = trim(adjustl(string))//C_NULL_CHAR diff --git a/src/vinit.F90 b/src/vinit.F90 index 79b5f2ac..665f50c8 100644 --- a/src/vinit.F90 +++ b/src/vinit.F90 @@ -54,7 +54,7 @@ subroutine vinit(temp, mass, vx, vy, vz) pom = 1 do iat = 1, natom ! Variance of distribution - sigma = dsqrt(TEMP / MASS(iat)) + sigma = sqrt(TEMP / MASS(iat)) ! Velocity vector of iat-th atom vx(iat, iw) = sigma * rans(pom) @@ -103,7 +103,7 @@ subroutine scale_velocities(vx, vy, vz) if (scaleveloc == 1 .and. temp_mom > 0.1E-10) then write (stdout, *) 'Scaling velocities to correct temperature.' - scal = dsqrt(temp / temp_mom) + scal = sqrt(temp / temp_mom) vx = vx * scal vy = vy * scal vz = vz * scal @@ -202,7 +202,7 @@ subroutine remove_rotations(x, y, z, vx, vy, vz, masses) Omx = Iinv(0) * Lx + Iinv(1) * Ly + Iinv(2) * Lz Omy = Iinv(3) * Lx + Iinv(4) * Ly + Iinv(5) * Lz Omz = Iinv(6) * Lx + Iinv(7) * Ly + Iinv(8) * Lz - Om_tot = dsqrt(Omx**2 + Omy**2 + Omz**2) + Om_tot = sqrt(Omx**2 + Omy**2 + Omz**2) ! We probably need to diagonalize I to get kinetic energy... !write(*,*)"Angular velocity:", Omx, Omy, Omz, Om_tot @@ -315,7 +315,7 @@ subroutine calc_angular_momentum(x, y, z, vx, vy, vz, masses, iw, Lx, Ly, Lz, L_ Lz = Lz + (x(iat, iw) * vy(iat, iw) - y(iat, iw) * vx(iat, iw)) * masses(iat) end do - L_tot = dsqrt(Lx**2 + Ly**2 + Lz**2) + L_tot = sqrt(Lx**2 + Ly**2 + Lz**2) end subroutine calc_angular_momentum subroutine constrainP(px, py, pz, constrained_atoms)