diff --git a/src/m_common/Thomas_Fermi_model.f90 b/src/m_common/Thomas_Fermi_model.f90
index 54c435e..24bb7af 100644
--- a/src/m_common/Thomas_Fermi_model.f90
+++ b/src/m_common/Thomas_Fermi_model.f90
@@ -38,10 +38,10 @@ module Thomas_Fermi_model
function Thomas_Fermi_core_electron_number(atomic_number, r_c) result(num_core_electron)
- real(8), intent(in) :: atomic_number, r_c
- real(8) num_core_electron, num_val_electron
- real(8) r_TF, x_TF, b, pi, rho_r, p_r
- real(8) :: r = 0.0d0, r_min = 1.0d-3, r_max = 25.0d0, dr = 1.0d-3
+ real*8, intent(in) :: atomic_number, r_c
+ real*8 num_core_electron, num_val_electron
+ real*8 r_TF, x_TF, b, pi, rho_r, p_r
+ real*8 :: r = 0.0d0, r_min = 1.0d-3, r_max = 25.0d0, dr = 1.0d-3
! define constant
pi = 3.1415926535d0
@@ -78,8 +78,8 @@ end function Thomas_Fermi_core_electron_number
function Gross_Dreizler(x) result(f)
- real(8), intent(in) :: x
- real(8) f
+ real*8, intent(in) :: x
+ real*8 f
f = 1.0d0/(1.0d0 + 1.4712d0*x - 0.4973d0*x**(3.0d0/2.0d0) + 0.3875d0*x**(2.0d0) + 0.002102d0*x**(3.0d0))
diff --git a/src/m_common/cell.f90 b/src/m_common/cell.f90
index 7e5ee1a..f190283 100644
--- a/src/m_common/cell.f90
+++ b/src/m_common/cell.f90
@@ -399,8 +399,11 @@ subroutine ApplyPBC(s, howmany)
end subroutine ApplyPBC
function map(x, cell_period)
- real(8) :: x, cell_period
- real(8) :: map
+ implicit none
+ ! argument variables
+ real*8, intent(in) :: x, cell_period
+ ! local variables
+ real*8 :: map
if (cell_period .eq. 0.d0) then
map = x
else
@@ -409,6 +412,7 @@ function map(x, cell_period)
end function map
function dmap(x, cell_period)
+ implicit none
real(8) :: x, cell_period
real(8) :: dmap
! dmap=dcos(x/cell_period)
@@ -420,6 +424,7 @@ function dmap(x, cell_period)
end function dmap
function ddmap(x, cell_period)
+ implicit none
real(8) :: x, cell_period
real(8) :: ddmap
if (cell_period .eq. 0.d0) then
@@ -430,7 +435,11 @@ function ddmap(x, cell_period)
end function ddmap
function map0(x)
- real(8) :: x, xc, map0
+ implicit none
+ ! argument variables
+ real*8, intent(in) :: x
+ ! local variables
+ real*8 xc, map0
integer p
! this function depend only on x and is such that f'=1 and f(1/2)=0
select case (case_map)
diff --git a/src/m_common/rotate_tools.f90 b/src/m_common/rotate_tools.f90
index 11072f8..9b01742 100644
--- a/src/m_common/rotate_tools.f90
+++ b/src/m_common/rotate_tools.f90
@@ -13,15 +13,13 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
-subroutine ruota_xyz(alpha, xrot, yrot, zrot, &
- & nion_1, rion_1, rion_2)
+subroutine ruota_xyz(alpha, xrot, yrot, zrot, nion_1, rion_1, rion_2)
implicit none
- real(8) u(3, 3)
- real(8), dimension(:, :), allocatable :: emme
+ real*8 u(3, 3)
+ real*8, dimension(:, :), allocatable :: emme
integer nion_1, i, j, ii
- real(8) alpha, xrot, yrot, zrot, &
- & rion_1(3, nion_1), rion_2(3, nion_1)
+ real*8 alpha, xrot, yrot, zrot, rion_1(3, nion_1), rion_2(3, nion_1)
call dscal(9, 0.d0, u(1, 1), 1)
call make_u(alpha, xrot, yrot, zrot, u)
@@ -43,19 +41,20 @@ subroutine ruota_xyz(alpha, xrot, yrot, zrot, &
end subroutine ruota_xyz
-subroutine ruota_molec(ipc, alpha, xrot, yrot, zrot, iesupr_1&
- &, dupr_1, ioptorb_1, nparam_1, nshell_1, ioptorb, nshell, nelorb, dupr_2)
+subroutine ruota_molec(ipc, alpha, xrot, yrot, zrot, iesupr_1, &
+ & dupr_1, ioptorb_1, nparam_1, nshell_1, &
+ & ioptorb, nshell, nelorb, dupr_2)
use constants, only: zzero, zone
implicit none
integer nshell_1, indpar, nelorb, nshell, shift, i, j, ii, iesupr_1, ipc
integer ioptorb(nshell), ioptorb_1(nshell_1), nparam_1(nshell_1)
- real(8) alpha, xrot, yrot, zrot
- real(8) dupr_1(ipc*iesupr_1)
- real(8) dupr_2(ipc*iesupr_1)
+ real*8 alpha, xrot, yrot, zrot
+ real*8 dupr_1(ipc*iesupr_1)
+ real*8 dupr_2(ipc*iesupr_1)
- real(8) u(3, 3), emmed(5, 5), emmef(7, 7), emmeg(9, 9)
- real(8), dimension(:, :), allocatable :: emme
+ real*8 u(3, 3), emmed(5, 5), emmef(7, 7), emmeg(9, 9)
+ real*8, dimension(:, :), allocatable :: emme
allocate (emme(nelorb, nelorb))
@@ -93,8 +92,9 @@ subroutine ruota_molec(ipc, alpha, xrot, yrot, zrot, iesupr_1&
end subroutine ruota_molec
subroutine ruota_lambda(ipc, ipf, alpha, xrot, yrot, zrot, &
- & ix_1, iy_1, detmat_1, nnozero_1, occ_1, nelcol, &
- & ioptorb_1, nshell_1, nnozero_2, ix_2, iy_2, detmat_2, symmagp)
+ & ix_1, iy_1, detmat_1, nnozero_1, occ_1, nelcol, &
+ & ioptorb_1, nshell_1, nnozero_2, ix_2, iy_2, detmat_2, &
+ & symmagp)
use allio, only: yes_hermite
use constants, only: zone, zzero
! INPUT
@@ -112,15 +112,15 @@ subroutine ruota_lambda(ipc, ipf, alpha, xrot, yrot, zrot, &
integer nshell_1, nelcol, neldiff, &
& ioptorb_1(nshell_1)
- real(8) alpha, xrot, yrot, zrot
+ real*8 alpha, xrot, yrot, zrot
integer ix_1(nnozero_1), iy_1(nnozero_1)
- real(8) detmat_1(*)
+ real*8 detmat_1(*)
integer ix_2(*), iy_2(*)
- real(8) detmat_2(*)
+ real*8 detmat_2(*)
integer i, j, ii
- real(8) u(3, 3), emmed(5, 5), emmef(7, 7), emmeg(9, 9)
- real(8), dimension(:, :), allocatable :: emme, lwork, lambda
+ real*8 u(3, 3), emmed(5, 5), emmef(7, 7), emmeg(9, 9)
+ real*8, dimension(:, :), allocatable :: emme, lwork, lambda
logical symmagp
allocate (lambda(ipc*occ_1, nelcol), &
diff --git a/src/m_common/save_jall_fn.f90 b/src/m_common/save_jall_fn.f90
index f8cea89..29f6c21 100644
--- a/src/m_common/save_jall_fn.f90
+++ b/src/m_common/save_jall_fn.f90
@@ -25,11 +25,12 @@ subroutine save_jall(yesfn, jastrowall_ee, winvjbar, winvjbarsz, winvj, psip)
use allio, only: nel, nelup, neldo, indt, indt4j, nelorbj, nelorbjh, iessz
implicit none
- ! input
- real(8), intent(inout) :: jastrowall_ee(nelup + neldo, nelup + neldo, 0:indt4j), psip(nel, nel)
- real(8), intent(in) :: winvj(max(nelorbjh, 1), 0:indt4j, *), winvjbar(*), winvjbarsz(*)
+ ! argument variables
+ real*8, intent(inout) :: jastrowall_ee(nelup + neldo, nelup + neldo, 0:indt4j), psip(nel, nel)
+ real*8, intent(in) :: winvj(max(nelorbjh, 1), 0:indt4j, *), winvjbar(*), winvjbarsz(*)
+
+ ! local variables
integer :: nelused
- ! local
integer :: i, k, nelorbj5
logical yesfn
diff --git a/src/m_common/scalevect.f90 b/src/m_common/scalevect.f90
index d50f133..2549997 100644
--- a/src/m_common/scalevect.f90
+++ b/src/m_common/scalevect.f90
@@ -16,8 +16,11 @@
subroutine scalevect(n, cellfat, vect)
use cell, only: map
implicit none
- integer n, i
- real*8 cellfat(3), vect(3*n)
+ integer, intent(in) :: n
+ real*8, intent(in) :: cellfat(3)
+ real*8, intent(inout) :: vect(3*n)
+
+ integer i
!$omp parallel do default(shared) private(i)
do i = 1, 3*n, 3
vect(i) = map(vect(i), cellfat(1))
diff --git a/src/m_common/slaterorb.f90 b/src/m_common/slaterorb.f90
index 27b0a42..d86239e 100644
--- a/src/m_common/slaterorb.f90
+++ b/src/m_common/slaterorb.f90
@@ -16,7 +16,7 @@
function slaterorb(ioptorb)
use constants, only: iflagerr
implicit none
- integer ioptorb
+ integer, intent(in) :: ioptorb
logical slaterorb
select case (ioptorb)
! case(34,10,12,28,57,80) ! Slater s
diff --git a/src/m_common/sub_comm.f90 b/src/m_common/sub_comm.f90
index 65dbe32..f96df56 100644
--- a/src/m_common/sub_comm.f90
+++ b/src/m_common/sub_comm.f90
@@ -20,14 +20,18 @@ module sub_comm
integer comm, parent
logical yesin
end type
+
contains
+
subroutine mpi_sub_comm_create(parent_comm, new_size, child, ierror)
implicit none
#ifdef PARALLEL
include 'mpif.h'
#endif
- integer parent_comm, new_size, ierror
- type(mpi_sub_comm) :: child
+ integer, intent(in) :: parent_comm
+ integer, intent(inout) :: new_size
+ integer, intent(inout) :: ierror ! intent(out)?
+ type(mpi_sub_comm), intent(inout) :: child
integer orig_group, sub_group, i
integer, dimension(:), allocatable :: ranks
@@ -74,8 +78,8 @@ end subroutine mpi_sub_comm_create
subroutine mpi_sub_comm_free(child, ierror)
implicit none
- integer ierror
- type(mpi_sub_comm) :: child
+ integer, intent(inout) :: ierror ! intent(out)?
+ type(mpi_sub_comm), intent(inout) :: child
#ifdef PARALLEL
include 'mpif.h'
diff --git a/src/m_common/symmetrize_agp.f90 b/src/m_common/symmetrize_agp.f90
index 7d12052..496d94a 100644
--- a/src/m_common/symmetrize_agp.f90
+++ b/src/m_common/symmetrize_agp.f90
@@ -18,15 +18,21 @@ subroutine symmetrizeagp(nnozero_c, nozero_c, jbradet, jbradetn, dsw&
&, symmagp, yes_hermite)
use constants, only: ipc, ipf, deps
use cell, only: cellscale, phase2pi
- use allio, only: rank, rion, yes_crystal, real_agp, molyes&
- &, pfaffup, kiontot
+ use allio, only: rank, rion, yes_crystal, real_agp, molyes, pfaffup, kiontot
implicit none
- integer nnozero_c, iessw0, iesswt, ii, jj, kk, ind, ix, iy&
- &, nelorb_c, nelcol_c, indc, iessw, nelorb_at, iyr, ierr, ndim
- integer nozero_c(*), itouch(*), jbradet(*), jbradetn(*)
- real*8 dsw(*), dsw0, detmat_c(ipc*nelorb_c*nelcol_c), scale_c(*)
- real*8 dswc(2), riondiff(3), cosphihalf, sinphihalf, cost, detr, deti
- logical symmagp, yes_hermite
+
+ ! argument variables
+ integer, intent(in) :: nnozero_c, iessw0, nelorb_c, nelcol_c, nelorb_at, &
+ & nozero_c(*), jbradet(*), jbradetn(*)
+ integer, intent(inout) :: itouch(*)
+ real*8, intent(inout) :: dsw(*)
+ real*8, intent(inout) :: detmat_c(ipc*nelorb_c*nelcol_c), scale_c(*)
+ logical, intent(in) :: symmagp, yes_hermite
+
+ ! local variables
+ integer ii, jj, kk, ind, indc, iessw, iesswt, ix, iy, iyr, ierr, ndim
+ real*8 dsw0, dswc(2), riondiff(3), cosphihalf, sinphihalf, cost, detr, deti
+
! Input detmat_c output detmat_c symmetrized
do ii = 1, nnozero_c
scale_c(ipc*(ii - 1) + 1:ipc*ii) = detmat_c(ipc*(nozero_c(ii) - 1) + 1:ipc*nozero_c(ii))
diff --git a/src/m_common/types.f90 b/src/m_common/types.f90
index 34200aa..1301362 100644
--- a/src/m_common/types.f90
+++ b/src/m_common/types.f90
@@ -35,13 +35,13 @@ module types
! ex. for Jastrow I have:
! vj,vju,winvj,jasmat,jasmat_sz,jasmat_c,jasmatsz_c
type wf_factor
- real(8), allocatable :: twobody_par(:)
- real(8), allocatable :: exps(:)
- real(8), allocatable :: bas_mat(:)
- real(8), allocatable :: exp_mat(:)
- real(8), allocatable :: exp_mat_sz(:)
- real(8), allocatable :: exp_mat_c(:)
- real(8), allocatable :: exp_mat_sz_c(:)
+ real*8, allocatable :: twobody_par(:)
+ real*8, allocatable :: exps(:)
+ real*8, allocatable :: bas_mat(:)
+ real*8, allocatable :: exp_mat(:)
+ real*8, allocatable :: exp_mat_sz(:)
+ real*8, allocatable :: exp_mat_c(:)
+ real*8, allocatable :: exp_mat_sz_c(:)
end type
end module types
diff --git a/src/m_common/update_jastrowall.f90 b/src/m_common/update_jastrowall.f90
index b8a0a6e..5f7877e 100644
--- a/src/m_common/update_jastrowall.f90
+++ b/src/m_common/update_jastrowall.f90
@@ -18,8 +18,14 @@
! -----------------------------------------
subroutine upjastrowall(nel, jastrowall, psip)
implicit none
- integer j, k, nel
- real*8 psip(nel, *), jastrowall(nel, *)
+
+ ! argument variables
+ integer, intent(in) :: nel
+ real*8, intent(in) :: psip(nel, *)
+ real*8, intent(inout) :: jastrowall(nel, *)
+
+ ! local variables
+ integer j, k
do j = 1, nel
do k = 1, nel
@@ -31,8 +37,14 @@ end subroutine upjastrowall
subroutine upjastrowall_sz(nel, nelup, jastrowall, psip)
implicit none
- integer j, k, nel, nelup
- real*8 psip(nel, *), jastrowall(nel, *)
+
+ ! argument variables
+ integer, intent(in) :: nel, nelup
+ real*8, intent(in) :: psip(nel, *)
+ real*8, intent(inout) :: jastrowall(nel, *)
+
+ ! local variables
+ integer j, k
do j = 1, nelup
do k = 1, nelup
@@ -59,8 +71,14 @@ end subroutine upjastrowall_sz
! -----------------------------------------
subroutine upjastrowallfat(nel, jastrowall, psipmu)
implicit none
- integer nel, i, j
- real*8 jastrowall(nel, *), psipmu(nel, *)
+
+ ! argument variables
+ integer, intent(in) :: nel
+ real*8, intent(in) :: psipmu(nel, *)
+ real*8, intent(inout) :: jastrowall(nel, *)
+
+ ! local variables
+ integer i, j
do i = 1, nel
do j = 1, nel
@@ -73,8 +91,14 @@ end subroutine upjastrowallfat
subroutine upjastrowallfat_sz(nel, nelup, jastrowall, psipmu)
implicit none
- integer nel, nelup, i, j
- real*8 jastrowall(nel, *), psipmu(nel, *)
+
+ ! argument variables
+ integer, intent(in) :: nel, nelup
+ real*8, intent(in) :: psipmu(nel, *)
+ real*8, intent(inout) :: jastrowall(nel, *)
+
+ ! local variables
+ integer i, j
do i = 1, nelup
do j = 1, nelup
@@ -103,8 +127,14 @@ end subroutine upjastrowallfat_sz
! ------------------------------------------
subroutine upjastrowallpsi(nel, jastrowall, psip)
implicit none
- integer j, k, nel
- real*8 psip(nel, *), jastrowall(nel, *)
+
+ ! argument variables
+ integer, intent(in) :: nel
+ real*8, intent(in) :: psip(nel, *)
+ real*8, intent(inout) :: jastrowall(nel, *)
+
+ ! local variables
+ integer j, k
do j = 1, nel
do k = 1, nel
@@ -117,8 +147,14 @@ end subroutine upjastrowallpsi
subroutine upjastrowallpsi_sz(nel, nelup, jastrowall, psip)
implicit none
- integer j, k, nel, nelup
- real*8 psip(nel, *), jastrowall(nel, *)
+
+ ! argument variables
+ integer, intent(in) :: nel, nelup
+ real*8, intent(in) :: psip(nel, *)
+ real*8, intent(inout) :: jastrowall(nel, *)
+
+ ! local variables
+ integer j, k
do j = 1, nelup
do k = j + 1, nelup
diff --git a/src/m_common/updatedwarp.f90 b/src/m_common/updatedwarp.f90
index 2a124d2..a8085b7 100644
--- a/src/m_common/updatedwarp.f90
+++ b/src/m_common/updatedwarp.f90
@@ -1,200 +1,204 @@
-!TL off
subroutine updatedwarp(nel, nion, kel, rion, rmu, r, kelelocb, kellogb&
&, rionelocb, rionlogb, cellb, celllb, force, pulay&
&, iespbc, warp, power, atom_number, warpmat, niong)
use allio, only: norm_metric
- use Cell
+ use Cell, only: CartesianToCrystal, map, cellscale, metric, dmap, car2cry
implicit none
- logical iespbc, warp
- integer nel, nion, i, j, k, n, niong, irefg, ireft
- real*8 kel(3, nel),rion(3, nion),rmu(3,max(nel,nion),nion), r(nel,nion)&
- &, kelelocb(3, nel), kellogb(3, nel), rionelocb(3, nion), rionlogb(3, nion)&
- &, force(3, nion), pulay(3, nion), sumw, sumdw, wfunc, wderiv, cellb(3), celllb(3)&
- &, derpul, wder, xmu(3), atom_number(*), warpmat(nion - niong, *)
- real*8 power, rc
- if(iespbc) then
+ ! argument variables
+ logical, intent(in) :: iespbc, warp
+ integer, intent(in) :: nel, nion, niong
+ real*8, intent(in) :: power, kel(3, nel), rion(3, nion), atom_number(*)
+ real*8, intent(inout) :: rmu(3, max(nel, nion), nion), r(nel, nion), &
+ &kelelocb(3, nel), kellogb(3, nel), rionelocb(3, nion), &
+ &rionlogb(3, nion), force(3, nion), pulay(3, nion), &
+ &cellb(3), celllb(3), warpmat(nion - niong, *)
+
+ ! local variables
+ integer i, j, k, n, irefg, ireft
+ real*8 sumw, sumdw, wfunc, wderiv, derpul, wder, xmu(3)
+
+ if (iespbc) then
do k = 1, nion
- if(atom_number(k).gt.0) then
+ if (atom_number(k) .gt. 0) then
do i = 1, nel
rmu(1, i, k) = kel(1, i) - rion(1, k)
rmu(2, i, k) = kel(2, i) - rion(2, k)
rmu(3, i, k) = kel(3, i) - rion(3, k)
- enddo
+ end do
call CartesianToCrystal(rmu(1, 1, k), nel)
do i = 1, nel
rmu(1, i, k) = map(rmu(1, i, k), cellscale(1))
rmu(2, i, k) = map(rmu(2, i, k), cellscale(2))
rmu(3, i, k) = map(rmu(3, i, k), cellscale(3))
- enddo
+ end do
else
do i = 1, nion
rmu(1, i, k) = rion(1, k) - rion(1, i)
rmu(2, i, k) = rion(2, k) - rion(2, i)
rmu(3, i, k) = rion(3, k) - rion(3, i)
- enddo
+ end do
call CartesianToCrystal(rmu(1, 1, k), nion)
do i = 1, nion
rmu(1, i, k) = map(rmu(1, i, k), cellscale(1))
rmu(2, i, k) = map(rmu(2, i, k), cellscale(2))
rmu(3, i, k) = map(rmu(3, i, k), cellscale(3))
- enddo
- endif
- enddo
+ end do
+ end if
+ end do
else
do k = 1, nion
- if(atom_number(k).gt.0) then
+ if (atom_number(k) .gt. 0) then
do i = 1, nel
rmu(1, i, k) = kel(1, i) - rion(1, k)
rmu(2, i, k) = kel(2, i) - rion(2, k)
rmu(3, i, k) = kel(3, i) - rion(3, k)
- enddo
+ end do
else
do i = 1, nion
rmu(1, i, k) = rion(1, k) - rion(1, i)
rmu(2, i, k) = rion(2, k) - rion(2, i)
rmu(3, i, k) = rion(3, k) - rion(3, i)
- enddo
- endif
- enddo
- endif
+ end do
+ end if
+ end do
+ end if
- if(iespbc) then
+ if (iespbc) then
do k = 1, nion
do i = 1, nel
r(i, k) = max(norm_metric(rmu(1, i, k), metric), 1d-9)
- enddo
- enddo
+ end do
+ end do
else
do k = 1, nion
do i = 1, nel
r(i, k) = max(dsqrt(sum(rmu(:, i, k)**2)), 1d-9)
- enddo
- enddo
- endif
+ end do
+ end do
+ end if
- if(iespbc) then
+ if (iespbc) then
do k = 1, nion
- if(atom_number(k).gt.0) then
+ if (atom_number(k) .gt. 0) then
do i = 1, nel
rmu(1, i, k) = kel(1, i) - rion(1, k)
rmu(2, i, k) = kel(2, i) - rion(2, k)
rmu(3, i, k) = kel(3, i) - rion(3, k)
- enddo
+ end do
call CartesianToCrystal(rmu(1, 1, k), nel)
do i = 1, nel
- xmu(1)=map(rmu(1,i,k),cellscale(1))
- xmu(2)=map(rmu(2,i,k),cellscale(2))
- xmu(3)=map(rmu(3,i,k),cellscale(3))
- rmu(1, i, k)=(metric(1,1)*xmu(1)+metric(1,2)*xmu(2)+metric(1,3)*xmu(3))&
- &*dmap(rmu(1, i, k), cellscale(1))
- rmu(2, i, k) =(metric(2,1)*xmu(1)+metric(2,2)*xmu(2)+metric(2,3)*xmu(3))&
- &*dmap(rmu(2, i, k), cellscale(2))
- rmu(3, i, k) =(metric(3,1)*xmu(1)+metric(3,2)*xmu(2)+metric(3,3)*xmu(3))&
- &*dmap(rmu(3, i, k), cellscale(3))
-! HERE rmu is the derivative of r vs r_cell times r
+ xmu(1) = map(rmu(1, i, k), cellscale(1))
+ xmu(2) = map(rmu(2, i, k), cellscale(2))
+ xmu(3) = map(rmu(3, i, k), cellscale(3))
+ rmu(1, i, k) = (metric(1, 1)*xmu(1) + metric(1, 2)*xmu(2) + metric(1, 3)*xmu(3)) &
+ *dmap(rmu(1, i, k), cellscale(1))
+ rmu(2, i, k) = (metric(2, 1)*xmu(1) + metric(2, 2)*xmu(2) + metric(2, 3)*xmu(3)) &
+ *dmap(rmu(2, i, k), cellscale(2))
+ rmu(3, i, k) = (metric(3, 1)*xmu(1) + metric(3, 2)*xmu(2) + metric(3, 3)*xmu(3)) &
+ *dmap(rmu(3, i, k), cellscale(3))
+! HERE rmu is the derivative of r vs r_cell times r
! chain rule for r_cell = car2cry x r_physical
! dr/dr_phisical = dr/dr_cell x dr_cell/dr_physical
- xmu(:)=rmu(:,i,k)
- rmu(1,i,k)=xmu(1)*car2cry(1,1)+xmu(2)*car2cry(2,1)+xmu(3)*car2cry(3,1)
- rmu(2,i,k)=xmu(1)*car2cry(1,2)+xmu(2)*car2cry(2,2)+xmu(3)*car2cry(3,2)
- rmu(3,i,k)=xmu(1)*car2cry(1,3)+xmu(2)*car2cry(2,3)+xmu(3)*car2cry(3,3)
- enddo
+ xmu(:) = rmu(:, i, k)
+ rmu(1, i, k) = xmu(1)*car2cry(1, 1) + xmu(2)*car2cry(2, 1) + xmu(3)*car2cry(3, 1)
+ rmu(2, i, k) = xmu(1)*car2cry(1, 2) + xmu(2)*car2cry(2, 2) + xmu(3)*car2cry(3, 2)
+ rmu(3, i, k) = xmu(1)*car2cry(1, 3) + xmu(2)*car2cry(2, 3) + xmu(3)*car2cry(3, 3)
+ end do
else
do i = 1, nion
rmu(1, i, k) = rion(1, k) - rion(1, i)
rmu(2, i, k) = rion(2, k) - rion(2, i)
rmu(3, i, k) = rion(3, k) - rion(3, i)
- enddo
+ end do
call CartesianToCrystal(rmu(1, 1, k), nion)
do i = 1, nion
- xmu(1)=map(rmu(1,i,k),cellscale(1))
- xmu(2)=map(rmu(2,i,k),cellscale(2))
- xmu(3)=map(rmu(3,i,k),cellscale(3))
- rmu(1, i, k)=(metric(1,1)*xmu(1)+metric(1,2)*xmu(2)+metric(1,3)*xmu(3))&
- &*dmap(rmu(1, i, k), cellscale(1))
- rmu(2, i, k) =(metric(2,1)*xmu(1)+metric(2,2)*xmu(2)+metric(2,3)*xmu(3))&
- &*dmap(rmu(2, i, k), cellscale(2))
- rmu(3, i, k) =(metric(3,1)*xmu(1)+metric(3,2)*xmu(2)+metric(3,3)*xmu(3))&
- &*dmap(rmu(3, i, k), cellscale(3))
-! HERE rmu is the derivative of r vs r_cell times r
+ xmu(1) = map(rmu(1, i, k), cellscale(1))
+ xmu(2) = map(rmu(2, i, k), cellscale(2))
+ xmu(3) = map(rmu(3, i, k), cellscale(3))
+ rmu(1, i, k) = (metric(1, 1)*xmu(1) + metric(1, 2)*xmu(2) + metric(1, 3)*xmu(3)) &
+ *dmap(rmu(1, i, k), cellscale(1))
+ rmu(2, i, k) = (metric(2, 1)*xmu(1) + metric(2, 2)*xmu(2) + metric(2, 3)*xmu(3)) &
+ *dmap(rmu(2, i, k), cellscale(2))
+ rmu(3, i, k) = (metric(3, 1)*xmu(1) + metric(3, 2)*xmu(2) + metric(3, 3)*xmu(3)) &
+ *dmap(rmu(3, i, k), cellscale(3))
+! HERE rmu is the derivative of r vs r_cell times r
! chain rule for r_cell = car2cry x r_physical
! dr/dr_phisical = dr/dr_cell x dr_cell/dr_physical
- xmu(:)=rmu(:,i,k)
- rmu(1,i,k)=xmu(1)*car2cry(1,1)+xmu(2)*car2cry(2,1)+xmu(3)*car2cry(3,1)
- rmu(2,i,k)=xmu(1)*car2cry(1,2)+xmu(2)*car2cry(2,2)+xmu(3)*car2cry(3,2)
- rmu(3,i,k)=xmu(1)*car2cry(1,3)+xmu(2)*car2cry(2,3)+xmu(3)*car2cry(3,3)
- enddo
- endif
- enddo
- endif
+ xmu(:) = rmu(:, i, k)
+ rmu(1, i, k) = xmu(1)*car2cry(1, 1) + xmu(2)*car2cry(2, 1) + xmu(3)*car2cry(3, 1)
+ rmu(2, i, k) = xmu(1)*car2cry(1, 2) + xmu(2)*car2cry(2, 2) + xmu(3)*car2cry(3, 2)
+ rmu(3, i, k) = xmu(1)*car2cry(1, 3) + xmu(2)*car2cry(2, 3) + xmu(3)*car2cry(3, 3)
+ end do
+ end if
+ end do
+ end if
! NB Obviously in the calculation of pressures all electron and ion
! positions have to be put in the same box with coordinates
! | r_i | < cellscale(i)/2, i=1,2,3.
-
force = rionelocb
pulay = rionlogb
- if(warp) then
+ if (warp) then
do n = 1, nel
sumw = 0.d0
- if(power.ne.0.d0) then
+ if (power .ne. 0.d0) then
do k = 1, nion
- if(atom_number(k).gt.0.d0) then
- wfunc = 1.d0 / r(n, k)**power
+ if (atom_number(k) .gt. 0.d0) then
+ wfunc = 1.d0/r(n, k)**power
sumw = sumw + wfunc
- endif
- enddo
- endif
+ end if
+ end do
+ end if
do j = 1, 3
sumdw = 0.d0
- if(power.ne.0.d0) then
+ if (power .ne. 0.d0) then
do k = 1, nion
! taken a factor 1/2 of the Jacobian into account
- if(atom_number(k).gt.0.d0) then
- wderiv = -power / 2.d0 / r(n, k)**(power + 2) * rmu(j, n, k)
+ if (atom_number(k) .gt. 0.d0) then
+ wderiv = -power/2.d0/r(n, k)**(power + 2)*rmu(j, n, k)
sumdw = sumdw + wderiv
- endif
- enddo
- endif
+ end if
+ end do
+ end if
do i = 1, nion
! taken a factor 1/2 of the Jacobian into account
! wderiv=-power/2.d0/r(n,i)**(power+2)*rmu(j,n,i)
! wfunc=1.d0/r(n,i)**power
- if(atom_number(i).gt.0.d0) then
+ if (atom_number(i) .gt. 0.d0) then
- if(power.eq.0.d0) then
+ if (power .eq. 0.d0) then
wder = 1.d0
derpul = 0.d0
else
- wderiv = -(power / 2.d0) / r(n, i)**(power + 2) * rmu(j, n, i)
- wfunc = 1.d0 / r(n, i)**power
+ wderiv = -(power/2.d0)/r(n, i)**(power + 2)*rmu(j, n, i)
+ wfunc = 1.d0/r(n, i)**power
! Jacobian contribution
- derpul = wderiv / sumw - wfunc / sumw**2 * sumdw
- wder = wfunc / sumw
- endif
+ derpul = wderiv/sumw - wfunc/sumw**2*sumdw
+ wder = wfunc/sumw
+ end if
! warp contribution
- pulay(j, i) = pulay(j, i) + wder * kellogb(j, n) + derpul
- force(j, i) = force(j, i) + wder * kelelocb(j, n)
- endif
- enddo
- enddo
- enddo
+ pulay(j, i) = pulay(j, i) + wder*kellogb(j, n) + derpul
+ force(j, i) = force(j, i) + wder*kelelocb(j, n)
+ end if
+ end do
+ end do
+ end do
irefg = 0
do n = 1, nion
- if(atom_number(n).le.0) then
+ if (atom_number(n) .le. 0) then
irefg = irefg + 1
sumw = 0.d0
- if(power.ne.0.d0) then
+ if (power .ne. 0.d0) then
do k = 1, nion
- if(atom_number(k).gt.0.d0) then
- wfunc = 1.d0 / r(k, n)**power
+ if (atom_number(k) .gt. 0.d0) then
+ wfunc = 1.d0/r(k, n)**power
sumw = sumw + wfunc
- endif
- enddo
- endif
+ end if
+ end do
+ end if
do j = 1, 3
ireft = 0
do i = 1, nion
@@ -202,57 +206,57 @@ subroutine updatedwarp(nel, nion, kel, rion, rmu, r, kelelocb, kellogb&
! wderiv=-power/2.d0/r(n,i)**(power+2)*rmu(j,n,i)
! wfunc=1.d0/r(n,i)**power
- if(atom_number(i).gt.0.d0) then
+ if (atom_number(i) .gt. 0.d0) then
ireft = ireft + 1
- if(power.ne.0.d0) then
- wfunc = 1.d0 / r(i, n)**power
- wder = wfunc / sumw
+ if (power .ne. 0.d0) then
+ wfunc = 1.d0/r(i, n)**power
+ wder = wfunc/sumw
else
wder = 1.d0
- endif
+ end if
! warp contribution
- pulay(j, i) = pulay(j, i) + wder * rionlogb(j, n)
- force(j, i) = force(j, i) + wder * rionelocb(j, n)
+ pulay(j, i) = pulay(j, i) + wder*rionlogb(j, n)
+ force(j, i) = force(j, i) + wder*rionelocb(j, n)
warpmat(ireft, irefg) = wder
- endif
- enddo
- enddo
- endif
- enddo
- endif ! endif warp
- if(iespbc) then
+ end if
+ end do
+ end do
+ end if
+ end do
+ end if ! endif warp
+ if (iespbc) then
! call ApplyPBC(rmu,nion*nel)
! press=cellb
! press_pulay=celllb
- rmu(:,1:nion,1)=rion(:,1:nion)
- call CartesianToCrystal(rmu,nion)
+ rmu(:, 1:nion, 1) = rion(:, 1:nion)
+ call CartesianToCrystal(rmu, nion)
! call ApplyPBC(rmu, nion)
do k = 1, nion
- cellb(1) = cellb(1) + rionelocb(1, k) * rmu(1, k,1) / cellscale(1)
- celllb(1) = celllb(1) + rionlogb(1, k) * rmu(1, k,1) / cellscale(1)
+ cellb(1) = cellb(1) + rionelocb(1, k)*rmu(1, k, 1)/cellscale(1)
+ celllb(1) = celllb(1) + rionlogb(1, k)*rmu(1, k, 1)/cellscale(1)
- cellb(2) = cellb(2) + rionelocb(2, k) * rmu(2, k,1) / cellscale(2)
- celllb(2) = celllb(2) + rionlogb(2, k) * rmu(2, k,1) / cellscale(2)
+ cellb(2) = cellb(2) + rionelocb(2, k)*rmu(2, k, 1)/cellscale(2)
+ celllb(2) = celllb(2) + rionlogb(2, k)*rmu(2, k, 1)/cellscale(2)
- cellb(3) = cellb(3) + rionelocb(3, k) * rmu(3, k,1) / cellscale(3)
- celllb(3) = celllb(3) + rionlogb(3, k) * rmu(3, k,1) / cellscale(3)
- enddo
- rmu(:,1:nel,1)=kel(:,1:nel)
+ cellb(3) = cellb(3) + rionelocb(3, k)*rmu(3, k, 1)/cellscale(3)
+ celllb(3) = celllb(3) + rionlogb(3, k)*rmu(3, k, 1)/cellscale(3)
+ end do
+ rmu(:, 1:nel, 1) = kel(:, 1:nel)
! call ApplyPBC(rmu, nel)
- call CartesianToCrystal(rmu,nel)
+ call CartesianToCrystal(rmu, nel)
do k = 1, nel
- cellb(1) = cellb(1) + kelelocb(1, k) * rmu(1, k,1) / cellscale(1)
- celllb(1) = celllb(1) + kellogb(1, k) * rmu(1, k,1) / cellscale(1)
+ cellb(1) = cellb(1) + kelelocb(1, k)*rmu(1, k, 1)/cellscale(1)
+ celllb(1) = celllb(1) + kellogb(1, k)*rmu(1, k, 1)/cellscale(1)
- cellb(2) = cellb(2) + kelelocb(2, k) * rmu(2, k,1) / cellscale(2)
- celllb(2) = celllb(2) + kellogb(2, k) * rmu(2, k,1) / cellscale(2)
+ cellb(2) = cellb(2) + kelelocb(2, k)*rmu(2, k, 1)/cellscale(2)
+ celllb(2) = celllb(2) + kellogb(2, k)*rmu(2, k, 1)/cellscale(2)
- cellb(3) = cellb(3) + kelelocb(3, k) * rmu(3, k,1) / cellscale(3)
- celllb(3) = celllb(3) + kellogb(3, k) * rmu(3, k,1) / cellscale(3)
- enddo
- endif
+ cellb(3) = cellb(3) + kelelocb(3, k)*rmu(3, k, 1)/cellscale(3)
+ celllb(3) = celllb(3) + kellogb(3, k)*rmu(3, k, 1)/cellscale(3)
+ end do
+ end if
return
end
diff --git a/src/m_common/upsim.f90 b/src/m_common/upsim.f90
index 270af11..4b24117 100644
--- a/src/m_common/upsim.f90
+++ b/src/m_common/upsim.f90
@@ -17,9 +17,15 @@ subroutine upsim(amat, ndim, ind, value, symmagp, ipf)
use allio, only: nelorb_at, pfaffup, kiontot
implicit none
- integer ndim, ndimh, ind, i, j, ipf
- real(8) amat(ndim, *), value
- logical symmagp
+ ! argument parameters
+ integer, intent(in) :: ndim, ind, ipf
+ real*8, intent(in) :: value
+ real*8, intent(inout) :: amat(ndim, *)
+ logical, intent(in) :: symmagp
+
+ ! local variables
+ integer :: i, j, ndimh
+
if (ipf .eq. 1) then
j = (ind - 1)/ndim + 1
i = ind - ndim*(j - 1)
@@ -59,9 +65,16 @@ end subroutine upsim
subroutine upsim_complex(amat, ndim, ind, value, symmagp, ipf)
use allio, only: yes_hermite, nelorb_at, pfaffup, kiontot
implicit none
- integer ndim, ndimh, ind, i, j, ipf
- complex*16 amat(ndim, *), value(*)
- logical symmagp
+
+ ! argument parameters
+ integer, intent(in) :: ndim, ind, ipf
+ complex*16, intent(in) :: value(*)
+ complex*16, intent(inout) :: amat(ndim, *)
+ logical, intent(in) :: symmagp
+
+ ! local variables
+ integer i, j, ndimh
+
! write(6,*) value
if (ipf .eq. 1) then
j = (ind - 1)/ndim + 1
diff --git a/src/m_common/upsimp.f90 b/src/m_common/upsimp.f90
index 17b30a5..5fc65a9 100644
--- a/src/m_common/upsimp.f90
+++ b/src/m_common/upsimp.f90
@@ -15,10 +15,18 @@
subroutine upsimp(amat, ndim, ind, value, symmagp, ipf)
use allio, only: nelorb_at, pfaffup, kiontot
+
implicit none
- integer ndim, ndimh, ind, i, j, ipf
- real*8 amat(ndim, *), value
- logical symmagp
+
+ ! argument parameters
+ integer, intent(in) :: ndim, ind, ipf
+ real*8, intent(in) :: value
+ real*8, intent(inout) :: amat(ndim, *)
+ logical, intent(in) :: symmagp
+
+ ! local variables
+ integer i, j, ndimh
+
if (ipf .eq. 1) then
j = (ind - 1)/ndim + 1
i = ind - ndim*(j - 1)
@@ -45,12 +53,21 @@ subroutine upsimp(amat, ndim, ind, value, symmagp, ipf)
end if
return
end
+
subroutine upsimp_complex(amat, ndim, ind, value, symmagp, ipf)
use allio, only: yes_hermite, nelorb_at, pfaffup, kiontot
+
implicit none
- integer ndim, ndimh, ind, ipf, i, j
- complex*16 amat(ndim, *), value(*)
- logical symmagp
+
+ ! argument parameters
+ integer, intent(in) :: ndim, ind, ipf
+ complex*16, intent(in) :: value(*)
+ complex*16, intent(inout) :: amat(ndim, *)
+ logical, intent(in) :: symmagp
+
+ ! local variables
+ integer i, j, ndimh
+
if (ipf .eq. 1) then
j = (ind - 1)/ndim + 1
i = ind - ndim*(j - 1)
diff --git a/src/m_common/upvpot_ei.f90 b/src/m_common/upvpot_ei.f90
index 259f12b..9b09698 100644
--- a/src/m_common/upvpot_ei.f90
+++ b/src/m_common/upvpot_ei.f90
@@ -15,9 +15,17 @@
subroutine upvpot_ei(rkel, zeta, rion, vpot, nion, LBox, epsvpot)
implicit none
- integer nel, k, nion
- real(8) rkel(3)
- real(8) zeta(nion), vpot, ngivej, LBox, rion(3, *), cost, epsvpot
+
+ ! argument parameters
+ integer nion
+ real*8, intent(in) :: rkel(3)
+ real*8, intent(in) :: zeta(nion), LBox, rion(3, *), epsvpot
+ real*8, intent(out) :: vpot
+
+ ! local variables
+ integer nel, k
+ real*8 cost, ngivej
+
if (epsvpot .eq. 0) then
vpot = 1.d0
else
diff --git a/src/m_common/upvpotaa.f90 b/src/m_common/upvpotaa.f90
index 66cdae0..60fcbfb 100644
--- a/src/m_common/upvpotaa.f90
+++ b/src/m_common/upvpotaa.f90
@@ -14,12 +14,21 @@
! along with this program. If not, see .
function upvpotaa(zeta, iond, nion, kappa, LBox)
- use dielectric
+ use dielectric, only: veps, rep_erfc, epsilon0
use allio, only: iond_cart, x_neigh, neigh, dist_shift, rank
+
implicit none
- integer nel, i, j, ii, jj, nion
- real(8) :: derfc, x_shift(3), cost_z
- real(8) zeta(nion), iond(nion, nion), kappa, LBox, pot_aa, upvpotaa, eself1b
+
+ ! argument parameters
+ integer, intent(in) :: nion
+ real*8, intent(in) :: zeta(nion), iond(nion, nion), kappa, LBox
+
+ ! local variables
+ integer nel, i, j, ii, jj
+ real*8 :: derfc, x_shift(3), cost_z
+ real*8 pot_aa, eself1b
+ real*8 upvpotaa
+
double precision, parameter :: PI = 3.14159265358979323846d0
pot_aa = 0.d0
@@ -28,7 +37,7 @@ function upvpotaa(zeta, iond, nion, kappa, LBox)
do i = 1, nion
do j = i + 1, nion
if (zeta(i)*zeta(j) .ne. 0.d0) &
- &pot_aa = pot_aa + 2.d0*zeta(i)*zeta(j)*veps(iond(i, j))
+ pot_aa = pot_aa + 2.d0*zeta(i)*zeta(j)*veps(iond(i, j))
end do
end do
! ******** EWALD REAL PART *************
@@ -45,10 +54,9 @@ function upvpotaa(zeta, iond, nion, kappa, LBox)
!$omp simd
#endif
do ii = 1, neigh
- dist_shift(ii) = dsqrt((iond_cart(1, jj)&
- & + x_neigh(ii, 1))**2 + (iond_cart(2, jj)&
- & + x_neigh(ii, 2))**2 + (iond_cart(3, jj)&
- & + x_neigh(ii, 3))**2)
+ dist_shift(ii) = dsqrt((iond_cart(1, jj) + x_neigh(ii, 1))**2 + &
+ (iond_cart(2, jj) + x_neigh(ii, 2))**2 + &
+ (iond_cart(3, jj) + x_neigh(ii, 3))**2)
dist_shift(ii) = cost_z*rep_erfc(dist_shift(ii), kappa)
end do
! pot_aa = pot_aa + cost_z*rep_erfc(dist_shift(ii),kappa)