PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 20:34:58 +01:00
Code Releases Activity

Fixed type error in localization

Browse Source
This commit is contained in:
Anthony Scemama 2022-11-19 16:41:44 +01:00
parent 5da97f2fbb
commit 50057251a7
5 changed files with 334 additions and 333 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
ocaml/qp_run.ml
View File

@ -110,7 +110,7 @@ let run slave ?prefix exe ezfio_file =
let task_thread = let task_thread =
let thread = let thread =
Thread.create ( fun () -> Thread.create ( fun () ->
TaskServer.run port_number ) TaskServer.run ~port:port_number )
in in
thread (); thread ();
in in

2
src/cis/cis.irp.f
View File

@ -79,6 +79,6 @@ subroutine run
call ezfio_set_cis_energy(CI_energy) call ezfio_set_cis_energy(CI_energy)
psi_coef = ci_eigenvectors psi_coef = ci_eigenvectors
SOFT_TOUCH psi_coef SOFT_TOUCH psi_coef
call save_wavefunction_truncated(thresh_save_wf) call save_wavefunction_truncated(save_threshold)
end end

7
src/determinants/EZFIO.cfg
View File

@ -136,9 +136,8 @@ doc: If |true|, discard any Slater determinants with an interaction smaller than
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
[save_threshold]
[thresh_save_wf]
type: Threshold type: Threshold
doc: Thresholds to save wave function doc: Cut-off to apply to the CI coefficients when the wave function is stored
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: 1.e-15 default: 1.e-14

385
src/mo_localization/localization.org
View File

File diff suppressed because it is too large Load Diff

271
src/mo_localization/localization_sub.irp.f
View File

@ -4,10 +4,10 @@
! Gradient: ! Gradient:
! qp_edit : ! qp_edit :
! | localization_method | method for the localization | ! | localization_method | method for the localization |
! Input: ! Input:
! | tmp_n | integer | Number of parameters in the MO subspace | ! | tmp_n | integer | Number of parameters in the MO subspace |
! | tmp_list_size | integer | Number of MOs in the mo_class we want to localize | ! | tmp_list_size | integer | Number of MOs in the mo_class we want to localize |
! | tmp_list(tmp_list_size) | integer | MOs in the mo_class | ! | tmp_list(tmp_list_size) | integer | MOs in the mo_class |
@ -20,7 +20,7 @@
subroutine gradient_localization(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) subroutine gradient_localization(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad)
include 'pi.h' include 'pi.h'
implicit none implicit none
@ -28,7 +28,7 @@ subroutine gradient_localization(tmp_n, tmp_list_size, tmp_list, v_grad, max_ele
BEGIN_DOC BEGIN_DOC
! Compute the gradient of the chosen localization method ! Compute the gradient of the chosen localization method
END_DOC END_DOC
integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad
@ -90,7 +90,7 @@ end
subroutine criterion_localization(tmp_list_size, tmp_list,criterion) subroutine criterion_localization(tmp_list_size, tmp_list,criterion)
include 'pi.h' include 'pi.h'
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -146,7 +146,7 @@ end
subroutine theta_localization(tmp_list, tmp_list_size, tmp_m_x, max_elem) subroutine theta_localization(tmp_list, tmp_list_size, tmp_m_x, max_elem)
include 'pi.h' include 'pi.h'
implicit none implicit none
@ -154,7 +154,7 @@ subroutine theta_localization(tmp_list, tmp_list_size, tmp_m_x, max_elem)
BEGIN_DOC BEGIN_DOC
! Compute the rotation angles between the MOs for the chosen localization method ! Compute the rotation angles between the MOs for the chosen localization method
END_DOC END_DOC
integer, intent(in) :: tmp_list_size, tmp_list(tmp_list_size) integer, intent(in) :: tmp_list_size, tmp_list(tmp_list_size)
double precision, intent(out) :: tmp_m_x(tmp_list_size,tmp_list_size), max_elem double precision, intent(out) :: tmp_m_x(tmp_list_size,tmp_list_size), max_elem
@ -170,7 +170,7 @@ subroutine theta_localization(tmp_list, tmp_list_size, tmp_m_x, max_elem)
end end
! Gradient ! Gradient
! Input: ! Input:
! | tmp_n | integer | Number of parameters in the MO subspace | ! | tmp_n | integer | Number of parameters in the MO subspace |
! | tmp_list_size | integer | Number of MOs in the mo_class we want to localize | ! | tmp_list_size | integer | Number of MOs in the mo_class we want to localize |
! | tmp_list(tmp_list_size) | integer | MOs in the mo_class | ! | tmp_list(tmp_list_size) | integer | MOs in the mo_class |
@ -188,13 +188,13 @@ end
subroutine gradient_FB(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) subroutine gradient_FB(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute the gradient for the Foster-Boys localization ! Compute the gradient for the Foster-Boys localization
END_DOC END_DOC
integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad
double precision, allocatable :: m_grad(:,:) double precision, allocatable :: m_grad(:,:)
@ -220,11 +220,11 @@ subroutine gradient_FB(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gr
+4d0 * mo_dipole_z(i,j) * (mo_dipole_z(i,i) - mo_dipole_z(j,j)) +4d0 * mo_dipole_z(i,j) * (mo_dipole_z(i,i) - mo_dipole_z(j,j))
enddo enddo
enddo enddo
! 2D -> 1D ! 2D -> 1D
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
call vec_to_mat_index(tmp_k,tmp_i,tmp_j) call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
v_grad(tmp_k) = m_grad(tmp_i,tmp_j) v_grad(tmp_k) = m_grad(tmp_i,tmp_j)
enddo enddo
! Maximum element in the gradient ! Maximum element in the gradient
@ -233,8 +233,8 @@ subroutine gradient_FB(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gr
if (ABS(v_grad(tmp_k)) > max_elem) then if (ABS(v_grad(tmp_k)) > max_elem) then
max_elem = ABS(v_grad(tmp_k)) max_elem = ABS(v_grad(tmp_k))
endif endif
enddo enddo
! Norm of the gradient ! Norm of the gradient
norm_grad = 0d0 norm_grad = 0d0
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
@ -243,7 +243,7 @@ subroutine gradient_FB(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gr
norm_grad = dsqrt(norm_grad) norm_grad = dsqrt(norm_grad)
print*, 'Maximal element in the gradient:', max_elem print*, 'Maximal element in the gradient:', max_elem
print*, 'Norm of the gradient:', norm_grad print*, 'Norm of the gradient:', norm_grad
! Deallocation ! Deallocation
deallocate(m_grad) deallocate(m_grad)
@ -261,7 +261,7 @@ end subroutine
! Gradient (OMP) ! Gradient (OMP)
subroutine gradient_FB_omp(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) subroutine gradient_FB_omp(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad)
use omp_lib use omp_lib
implicit none implicit none
@ -269,7 +269,7 @@ subroutine gradient_FB_omp(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, nor
BEGIN_DOC BEGIN_DOC
! Compute the gradient for the Foster-Boys localization ! Compute the gradient for the Foster-Boys localization
END_DOC END_DOC
integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad
double precision, allocatable :: m_grad(:,:) double precision, allocatable :: m_grad(:,:)
@ -310,8 +310,8 @@ subroutine gradient_FB_omp(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, nor
!$OMP DO !$OMP DO
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
call vec_to_mat_index(tmp_k,tmp_i,tmp_j) call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
v_grad(tmp_k) = m_grad(tmp_i,tmp_j) v_grad(tmp_k) = m_grad(tmp_i,tmp_j)
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
@ -324,7 +324,7 @@ subroutine gradient_FB_omp(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, nor
if (ABS(v_grad(tmp_k)) > max_elem) then if (ABS(v_grad(tmp_k)) > max_elem) then
max_elem = ABS(v_grad(tmp_k)) max_elem = ABS(v_grad(tmp_k))
endif endif
enddo enddo
! Norm of the gradient ! Norm of the gradient
norm_grad = 0d0 norm_grad = 0d0
@ -334,7 +334,7 @@ subroutine gradient_FB_omp(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, nor
norm_grad = dsqrt(norm_grad) norm_grad = dsqrt(norm_grad)
print*, 'Maximal element in the gradient:', max_elem print*, 'Maximal element in the gradient:', max_elem
print*, 'Norm of the gradient:', norm_grad print*, 'Norm of the gradient:', norm_grad
! Deallocation ! Deallocation
deallocate(m_grad) deallocate(m_grad)
@ -349,7 +349,7 @@ subroutine gradient_FB_omp(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, nor
end subroutine end subroutine
! Hessian ! Hessian
! Output: ! Output:
! | H(tmp_n,tmp_n) | double precision | Gradient in the subspace | ! | H(tmp_n,tmp_n) | double precision | Gradient in the subspace |
@ -367,7 +367,7 @@ end subroutine
subroutine hessian_FB(tmp_n, tmp_list_size, tmp_list, H) subroutine hessian_FB(tmp_n, tmp_list_size, tmp_list, H)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute the diagonal hessian for the Foster-Boys localization ! Compute the diagonal hessian for the Foster-Boys localization
END_DOC END_DOC
@ -377,7 +377,7 @@ subroutine hessian_FB(tmp_n, tmp_list_size, tmp_list, H)
double precision, allocatable :: beta(:,:) double precision, allocatable :: beta(:,:)
integer :: i,j,tmp_k,tmp_i, tmp_j integer :: i,j,tmp_k,tmp_i, tmp_j
double precision :: max_elem, t1,t2,t3 double precision :: max_elem, t1,t2,t3
print*,'' print*,''
print*,'---hessian_FB---' print*,'---hessian_FB---'
print*,'' print*,''
@ -387,7 +387,7 @@ subroutine hessian_FB(tmp_n, tmp_list_size, tmp_list, H)
! Allocation ! Allocation
allocate(beta(tmp_list_size,tmp_list_size)) allocate(beta(tmp_list_size,tmp_list_size))
! Calculation ! Calculation
do tmp_j = 1, tmp_list_size do tmp_j = 1, tmp_list_size
j = tmp_list(tmp_j) j = tmp_list(tmp_j)
@ -405,7 +405,7 @@ subroutine hessian_FB(tmp_n, tmp_list_size, tmp_list, H)
call vec_to_mat_index(tmp_k,tmp_i,tmp_j) call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
H(tmp_k,tmp_k) = 4d0 * beta(tmp_i, tmp_j) H(tmp_k,tmp_k) = 4d0 * beta(tmp_i, tmp_j)
enddo enddo
! Min elem ! Min elem
max_elem = 0d0 max_elem = 0d0
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
@ -423,7 +423,7 @@ subroutine hessian_FB(tmp_n, tmp_list_size, tmp_list, H)
endif endif
enddo enddo
print*, 'Max elem H:', max_elem print*, 'Max elem H:', max_elem
! Near 0 ! Near 0
max_elem = 1d10 max_elem = 1d10
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
@ -435,7 +435,7 @@ subroutine hessian_FB(tmp_n, tmp_list_size, tmp_list, H)
! Deallocation ! Deallocation
deallocate(beta) deallocate(beta)
call wall_time(t2) call wall_time(t2)
t3 = t2 - t1 t3 = t2 - t1
print*,'Time in hessian_FB:', t3 print*,'Time in hessian_FB:', t3
@ -461,7 +461,7 @@ subroutine hessian_FB_omp(tmp_n, tmp_list_size, tmp_list, H)
double precision, allocatable :: beta(:,:) double precision, allocatable :: beta(:,:)
integer :: i,j,tmp_k,tmp_i,tmp_j integer :: i,j,tmp_k,tmp_i,tmp_j
double precision :: max_elem, t1,t2,t3 double precision :: max_elem, t1,t2,t3
print*,'' print*,''
print*,'---hessian_FB_omp---' print*,'---hessian_FB_omp---'
print*,'' print*,''
@ -479,7 +479,7 @@ subroutine hessian_FB_omp(tmp_n, tmp_list_size, tmp_list, H)
!$OMP SHARED(tmp_n,tmp_list_size,beta,H,mo_dipole_x,mo_dipole_y,mo_dipole_z,tmp_list) & !$OMP SHARED(tmp_n,tmp_list_size,beta,H,mo_dipole_x,mo_dipole_y,mo_dipole_z,tmp_list) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
! Calculation ! Calculation
!$OMP DO !$OMP DO
do tmp_j = 1, tmp_list_size do tmp_j = 1, tmp_list_size
@ -497,11 +497,11 @@ subroutine hessian_FB_omp(tmp_n, tmp_list_size, tmp_list, H)
!$OMP DO !$OMP DO
do j = 1, tmp_n do j = 1, tmp_n
do i = 1, tmp_n do i = 1, tmp_n
H(i,j) = 0d0 H(i,j) = 0d0
enddo enddo
enddo enddo
!$OMP END DO !$OMP END DO
! Diagonalm of the hessian ! Diagonalm of the hessian
!$OMP DO !$OMP DO
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
@ -509,7 +509,7 @@ subroutine hessian_FB_omp(tmp_n, tmp_list_size, tmp_list, H)
H(tmp_k,tmp_k) = 4d0 * beta(tmp_i, tmp_j) H(tmp_k,tmp_k) = 4d0 * beta(tmp_i, tmp_j)
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
call omp_set_max_active_levels(4) call omp_set_max_active_levels(4)
@ -531,7 +531,7 @@ subroutine hessian_FB_omp(tmp_n, tmp_list_size, tmp_list, H)
endif endif
enddo enddo
print*, 'Max elem H:', max_elem print*, 'Max elem H:', max_elem
! Near 0 ! Near 0
max_elem = 1d10 max_elem = 1d10
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
@ -543,7 +543,7 @@ subroutine hessian_FB_omp(tmp_n, tmp_list_size, tmp_list, H)
! Deallocation ! Deallocation
deallocate(beta) deallocate(beta)
call wall_time(t2) call wall_time(t2)
t3 = t2 - t1 t3 = t2 - t1
print*,'Time in hessian_FB_omp:', t3 print*,'Time in hessian_FB_omp:', t3
@ -567,12 +567,12 @@ subroutine grad_pipek(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gra
integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad
double precision, allocatable :: m_grad(:,:), tmp_int(:,:) double precision, allocatable :: m_grad(:,:), tmp_int(:,:)
integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho
! Allocation ! Allocation
allocate(m_grad(tmp_list_size, tmp_list_size), tmp_int(tmp_list_size, tmp_list_size)) allocate(m_grad(tmp_list_size, tmp_list_size), tmp_int(tmp_list_size, tmp_list_size))
! Initialization ! Initialization
m_grad = 0d0 m_grad = 0d0
do a = 1, nucl_num ! loop over the nuclei do a = 1, nucl_num ! loop over the nuclei
@ -580,12 +580,12 @@ subroutine grad_pipek(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gra
! Loop over the MOs of the a given mo_class to compute <i|P_a|j> ! Loop over the MOs of the a given mo_class to compute <i|P_a|j>
do tmp_j = 1, tmp_list_size do tmp_j = 1, tmp_list_size
j = tmp_list(tmp_j) j = tmp_list(tmp_j)
do tmp_i = 1, tmp_list_size do tmp_i = 1, tmp_list_size
i = tmp_list(tmp_i) i = tmp_list(tmp_i)
do rho = 1, ao_num ! loop over all the AOs do rho = 1, ao_num ! loop over all the AOs
do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
mu = nucl_aos(a,b) ! AO centered on atom a mu = nucl_aos(a,b) ! AO centered on atom a
tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) & tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j)) + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
@ -609,7 +609,7 @@ subroutine grad_pipek(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gra
! 2D -> 1D ! 2D -> 1D
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
call vec_to_mat_index(tmp_k,tmp_i,tmp_j) call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
v_grad(tmp_k) = m_grad(tmp_i,tmp_j) v_grad(tmp_k) = m_grad(tmp_i,tmp_j)
enddo enddo
! Maximum element in the gradient ! Maximum element in the gradient
@ -654,7 +654,7 @@ end subroutine grad_pipek
! $\sum_{\mu \in A}$ -> sum over the AOs which belongs to atom A ! $\sum_{\mu \in A}$ -> sum over the AOs which belongs to atom A
! $c^t$ -> expansion coefficient of orbital |t> ! $c^t$ -> expansion coefficient of orbital |t>
! Input: ! Input:
! | tmp_n | integer | Number of parameters in the MO subspace | ! | tmp_n | integer | Number of parameters in the MO subspace |
! | tmp_list_size | integer | Number of MOs in the mo_class we want to localize | ! | tmp_list_size | integer | Number of MOs in the mo_class we want to localize |
! | tmp_list(tmp_list_size) | integer | MOs in the mo_class | ! | tmp_list(tmp_list_size) | integer | MOs in the mo_class |
@ -689,7 +689,7 @@ subroutine gradient_PM(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gr
BEGIN_DOC BEGIN_DOC
! Compute gradient for the Pipek-Mezey localization ! Compute gradient for the Pipek-Mezey localization
END_DOC END_DOC
integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad
double precision, allocatable :: m_grad(:,:), tmp_int(:,:), CS(:,:), tmp_mo_coef(:,:), tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:) double precision, allocatable :: m_grad(:,:), tmp_int(:,:), CS(:,:), tmp_mo_coef(:,:), tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:)
@ -713,12 +713,12 @@ subroutine gradient_PM(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gr
do j = 1, ao_num do j = 1, ao_num
tmp_mo_coef(j,tmp_i) = mo_coef(j,i) tmp_mo_coef(j,tmp_i) = mo_coef(j,i)
enddo enddo
enddo enddo
call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1)) call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
m_grad = 0d0 m_grad = 0d0
do a = 1, nucl_num ! loop over the nuclei do a = 1, nucl_num ! loop over the nuclei
@ -748,7 +748,7 @@ subroutine gradient_PM(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gr
enddo enddo
enddo enddo
do b = 1, nucl_n_aos(a) do b = 1, nucl_n_aos(a)
mu = nucl_aos(a,b) mu = nucl_aos(a,b)
do tmp_i = 1, tmp_list_size do tmp_i = 1, tmp_list_size
@ -756,14 +756,14 @@ subroutine gradient_PM(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gr
tmp_CS(tmp_i,b) = CS(tmp_i,mu) tmp_CS(tmp_i,b) = CS(tmp_i,mu)
enddo enddo
enddo enddo
call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1)) call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1))
do tmp_j = 1, tmp_list_size do tmp_j = 1, tmp_list_size
do tmp_i = 1, tmp_list_size do tmp_i = 1, tmp_list_size
tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i)) tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i))
enddo enddo
enddo enddo
@ -783,7 +783,7 @@ subroutine gradient_PM(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gr
! 2D -> 1D ! 2D -> 1D
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
call vec_to_mat_index(tmp_k,tmp_i,tmp_j) call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
v_grad(tmp_k) = m_grad(tmp_i,tmp_j) v_grad(tmp_k) = m_grad(tmp_i,tmp_j)
enddo enddo
! Maximum element in the gradient ! Maximum element in the gradient
@ -792,7 +792,7 @@ subroutine gradient_PM(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_gr
if (ABS(v_grad(tmp_k)) > max_elem) then if (ABS(v_grad(tmp_k)) > max_elem) then
max_elem = ABS(v_grad(tmp_k)) max_elem = ABS(v_grad(tmp_k))
endif endif
enddo enddo
! Norm of the gradient ! Norm of the gradient
norm_grad = 0d0 norm_grad = 0d0
@ -832,7 +832,7 @@ subroutine hess_pipek(tmp_n, tmp_list_size, tmp_list, H)
double precision, allocatable :: beta(:,:),tmp_int(:,:) double precision, allocatable :: beta(:,:),tmp_int(:,:)
integer :: i,j,tmp_k,tmp_i, tmp_j, a,b,rho,mu integer :: i,j,tmp_k,tmp_i, tmp_j, a,b,rho,mu
double precision :: max_elem double precision :: max_elem
! Allocation ! Allocation
allocate(beta(tmp_list_size,tmp_list_size),tmp_int(tmp_list_size,tmp_list_size)) allocate(beta(tmp_list_size,tmp_list_size),tmp_int(tmp_list_size,tmp_list_size))
@ -853,7 +853,7 @@ subroutine hess_pipek(tmp_n, tmp_list_size, tmp_list, H)
+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j)) + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
enddo enddo
enddo enddo
enddo enddo
enddo enddo
@ -865,7 +865,7 @@ subroutine hess_pipek(tmp_n, tmp_list_size, tmp_list, H)
enddo enddo
enddo enddo
enddo enddo
H = 0d0 H = 0d0
@ -873,7 +873,7 @@ subroutine hess_pipek(tmp_n, tmp_list_size, tmp_list, H)
call vec_to_mat_index(tmp_k,tmp_i,tmp_j) call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
H(tmp_k,tmp_k) = 4d0 * beta(tmp_i, tmp_j) H(tmp_k,tmp_k) = 4d0 * beta(tmp_i, tmp_j)
enddo enddo
! max_elem = 0d0 ! max_elem = 0d0
! do tmp_k = 1, tmp_n ! do tmp_k = 1, tmp_n
! if (H(tmp_k,tmp_k) < max_elem) then ! if (H(tmp_k,tmp_k) < max_elem) then
@ -889,7 +889,7 @@ subroutine hess_pipek(tmp_n, tmp_list_size, tmp_list, H)
! endif ! endif
! enddo ! enddo
! print*, 'Max elem H:', max_elem ! print*, 'Max elem H:', max_elem
! !
! max_elem = 1d10 ! max_elem = 1d10
! do tmp_k = 1, tmp_n ! do tmp_k = 1, tmp_n
! if (ABS(H(tmp_k,tmp_k)) < ABS(max_elem)) then ! if (ABS(H(tmp_k,tmp_k)) < ABS(max_elem)) then
@ -925,7 +925,7 @@ end
subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H) subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute diagonal hessian for the Pipek-Mezey localization ! Compute diagonal hessian for the Pipek-Mezey localization
END_DOC END_DOC
@ -935,7 +935,7 @@ subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
double precision, allocatable :: beta(:,:),tmp_int(:,:),CS(:,:),tmp_mo_coef(:,:),tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:) double precision, allocatable :: beta(:,:),tmp_int(:,:),CS(:,:),tmp_mo_coef(:,:),tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:)
integer :: i,j,tmp_k,tmp_i, tmp_j, a,b,rho,mu integer :: i,j,tmp_k,tmp_i, tmp_j, a,b,rho,mu
double precision :: max_elem, t1,t2,t3 double precision :: max_elem, t1,t2,t3
print*,'' print*,''
print*,'---hessian_PM---' print*,'---hessian_PM---'
print*,'' print*,''
@ -953,12 +953,12 @@ subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
do j = 1, ao_num do j = 1, ao_num
tmp_mo_coef(j,tmp_i) = mo_coef(j,i) tmp_mo_coef(j,tmp_i) = mo_coef(j,i)
enddo enddo
enddo enddo
call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1)) call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
do a = 1, nucl_num ! loop over the nuclei do a = 1, nucl_num ! loop over the nuclei
tmp_int = 0d0 tmp_int = 0d0
@ -975,7 +975,7 @@ subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
! enddo ! enddo
! enddo ! enddo
!enddo !enddo
allocate(tmp_mo_coef2(nucl_n_aos(a),tmp_list_size),tmp_CS(tmp_list_size,nucl_n_aos(a))) allocate(tmp_mo_coef2(nucl_n_aos(a),tmp_list_size),tmp_CS(tmp_list_size,nucl_n_aos(a)))
do tmp_i = 1, tmp_list_size do tmp_i = 1, tmp_list_size
@ -986,7 +986,7 @@ subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
enddo enddo
enddo enddo
do b = 1, nucl_n_aos(a) do b = 1, nucl_n_aos(a)
mu = nucl_aos(a,b) mu = nucl_aos(a,b)
do tmp_i = 1, tmp_list_size do tmp_i = 1, tmp_list_size
@ -994,14 +994,14 @@ subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
tmp_CS(tmp_i,b) = CS(tmp_i,mu) tmp_CS(tmp_i,b) = CS(tmp_i,mu)
enddo enddo
enddo enddo
call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1)) call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1))
do tmp_j = 1, tmp_list_size do tmp_j = 1, tmp_list_size
do tmp_i = 1, tmp_list_size do tmp_i = 1, tmp_list_size
tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i)) tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i))
enddo enddo
enddo enddo
@ -1016,7 +1016,7 @@ subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
enddo enddo
enddo enddo
enddo enddo
H = 0d0 H = 0d0
@ -1024,7 +1024,7 @@ subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
call vec_to_mat_index(tmp_k,tmp_i,tmp_j) call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
H(tmp_k,tmp_k) = 4d0 * beta(tmp_i, tmp_j) H(tmp_k,tmp_k) = 4d0 * beta(tmp_i, tmp_j)
enddo enddo
max_elem = 0d0 max_elem = 0d0
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
if (H(tmp_k,tmp_k) < max_elem) then if (H(tmp_k,tmp_k) < max_elem) then
@ -1040,7 +1040,7 @@ subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
endif endif
enddo enddo
print*, 'Max elem H:', max_elem print*, 'Max elem H:', max_elem
max_elem = 1d10 max_elem = 1d10
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
if (ABS(H(tmp_k,tmp_k)) < ABS(max_elem)) then if (ABS(H(tmp_k,tmp_k)) < ABS(max_elem)) then
@ -1067,18 +1067,18 @@ end
subroutine compute_crit_pipek(criterion) subroutine compute_crit_pipek(criterion)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute the Pipek-Mezey localization criterion ! Compute the Pipek-Mezey localization criterion
END_DOC END_DOC
double precision, intent(out) :: criterion double precision, intent(out) :: criterion
double precision, allocatable :: tmp_int(:,:) double precision, allocatable :: tmp_int(:,:)
integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho
! Allocation ! Allocation
allocate(tmp_int(mo_num, mo_num)) allocate(tmp_int(mo_num, mo_num))
criterion = 0d0 criterion = 0d0
do a = 1, nucl_num ! loop over the nuclei do a = 1, nucl_num ! loop over the nuclei
@ -1093,16 +1093,16 @@ subroutine compute_crit_pipek(criterion)
+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,i)) + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,i))
enddo enddo
enddo enddo
enddo enddo
do i = 1, mo_num do i = 1, mo_num
criterion = criterion + tmp_int(i,i)**2 criterion = criterion + tmp_int(i,i)**2
enddo enddo
enddo enddo
criterion = - criterion criterion = - criterion
deallocate(tmp_int) deallocate(tmp_int)
@ -1114,7 +1114,7 @@ end
! \begin{align*} ! \begin{align*}
! \mathcal{P} = \sum_{i=1}^n \sum_{A=1}^N \left[ <i|P_A|i> \right]^2 ! \mathcal{P} = \sum_{i=1}^n \sum_{A=1}^N \left[ <i|P_A|i> \right]^2
! \end{align*} ! \end{align*}
! with ! with
! \begin{align*} ! \begin{align*}
! <s|P_A|t> = \frac{1}{2} \sum_{\rho} \sum_{\mu \in A} \left[ c_{\rho}^{s*} S_{\rho \nu} c_{\mu}^{t} +c_{\mu}^{s*} S_{\mu \rho} c_{\rho}^t \right] ! <s|P_A|t> = \frac{1}{2} \sum_{\rho} \sum_{\mu \in A} \left[ c_{\rho}^{s*} S_{\rho \nu} c_{\mu}^{t} +c_{\mu}^{s*} S_{\mu \rho} c_{\rho}^t \right]
! \end{align*} ! \end{align*}
@ -1123,7 +1123,7 @@ end
subroutine criterion_PM(tmp_list_size,tmp_list,criterion) subroutine criterion_PM(tmp_list_size,tmp_list,criterion)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute the Pipek-Mezey localization criterion ! Compute the Pipek-Mezey localization criterion
END_DOC END_DOC
@ -1132,18 +1132,18 @@ subroutine criterion_PM(tmp_list_size,tmp_list,criterion)
double precision, intent(out) :: criterion double precision, intent(out) :: criterion
double precision, allocatable :: tmp_int(:,:),CS(:,:) double precision, allocatable :: tmp_int(:,:),CS(:,:)
integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho
print*,'' print*,''
print*,'---criterion_PM---' print*,'---criterion_PM---'
! Allocation ! Allocation
allocate(tmp_int(tmp_list_size, tmp_list_size),CS(mo_num,ao_num)) allocate(tmp_int(tmp_list_size, tmp_list_size),CS(mo_num,ao_num))
! Initialization ! Initialization
criterion = 0d0 criterion = 0d0
call dgemm('T','N',mo_num,ao_num,ao_num,1d0,mo_coef,size(mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1)) call dgemm('T','N',mo_num,ao_num,ao_num,1d0,mo_coef,size(mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
do a = 1, nucl_num ! loop over the nuclei do a = 1, nucl_num ! loop over the nuclei
tmp_int = 0d0 tmp_int = 0d0
@ -1151,7 +1151,7 @@ subroutine criterion_PM(tmp_list_size,tmp_list,criterion)
i = tmp_list(tmp_i) i = tmp_list(tmp_i)
do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
mu = nucl_aos(a,b) mu = nucl_aos(a,b)
tmp_int(tmp_i,tmp_i) = tmp_int(tmp_i,tmp_i) + 0.5d0 * (CS(i,mu) * mo_coef(mu,i) + mo_coef(mu,i) * CS(i,mu)) tmp_int(tmp_i,tmp_i) = tmp_int(tmp_i,tmp_i) + 0.5d0 * (CS(i,mu) * mo_coef(mu,i) + mo_coef(mu,i) * CS(i,mu))
! (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) & ! (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
@ -1165,8 +1165,8 @@ subroutine criterion_PM(tmp_list_size,tmp_list,criterion)
enddo enddo
enddo enddo
criterion = - criterion criterion = - criterion
deallocate(tmp_int,CS) deallocate(tmp_int,CS)
@ -1180,11 +1180,11 @@ end
subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion) subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute the Pipek-Mezey localization criterion ! Compute the Pipek-Mezey localization criterion
END_DOC END_DOC
integer, intent(in) :: tmp_list_size, tmp_list(tmp_list_size) integer, intent(in) :: tmp_list_size, tmp_list(tmp_list_size)
double precision, intent(out) :: criterion double precision, intent(out) :: criterion
double precision, allocatable :: tmp_int(:,:), CS(:,:), tmp_mo_coef(:,:), tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:) double precision, allocatable :: tmp_int(:,:), CS(:,:), tmp_mo_coef(:,:), tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:)
@ -1193,7 +1193,7 @@ subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion)
print*,'' print*,''
print*,'---criterion_PM_v3---' print*,'---criterion_PM_v3---'
call wall_time(t1) call wall_time(t1)
! Allocation ! Allocation
@ -1208,16 +1208,16 @@ subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion)
do j = 1, ao_num do j = 1, ao_num
tmp_mo_coef(j,tmp_i) = mo_coef(j,i) tmp_mo_coef(j,tmp_i) = mo_coef(j,i)
enddo enddo
enddo enddo
! ao_overlap(ao_num,ao_num) ! ao_overlap(ao_num,ao_num)
! mo_coef(ao_num,mo_num) ! mo_coef(ao_num,mo_num)
call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1)) call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
do a = 1, nucl_num ! loop over the nuclei do a = 1, nucl_num ! loop over the nuclei
do j = 1, tmp_list_size do j = 1, tmp_list_size
do i = 1, tmp_list_size do i = 1, tmp_list_size
tmp_int(i,j) = 0d0 tmp_int(i,j) = 0d0
@ -1239,7 +1239,7 @@ subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion)
!enddo !enddo
allocate(tmp_mo_coef2(nucl_n_aos(a),tmp_list_size),tmp_CS(tmp_list_size,nucl_n_aos(a))) allocate(tmp_mo_coef2(nucl_n_aos(a),tmp_list_size),tmp_CS(tmp_list_size,nucl_n_aos(a)))
do tmp_i = 1, tmp_list_size do tmp_i = 1, tmp_list_size
do b = 1, nucl_n_aos(a) do b = 1, nucl_n_aos(a)
mu = nucl_aos(a,b) mu = nucl_aos(a,b)
@ -1248,7 +1248,7 @@ subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion)
enddo enddo
enddo enddo
do b = 1, nucl_n_aos(a) do b = 1, nucl_n_aos(a)
mu = nucl_aos(a,b) mu = nucl_aos(a,b)
do tmp_i = 1, tmp_list_size do tmp_i = 1, tmp_list_size
@ -1256,15 +1256,15 @@ subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion)
tmp_CS(tmp_i,b) = CS(tmp_i,mu) tmp_CS(tmp_i,b) = CS(tmp_i,mu)
enddo enddo
enddo enddo
call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1)) call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1))
! Integrals ! Integrals
do tmp_j = 1, tmp_list_size do tmp_j = 1, tmp_list_size
do tmp_i = 1, tmp_list_size do tmp_i = 1, tmp_list_size
tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i)) tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i))
enddo enddo
enddo enddo
@ -1278,7 +1278,7 @@ subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion)
enddo enddo
criterion = - criterion criterion = - criterion
deallocate(tmp_int,CS,tmp_accu,tmp_mo_coef) deallocate(tmp_int,CS,tmp_accu,tmp_mo_coef)
@ -1477,14 +1477,14 @@ subroutine theta_FB(l, n, m_x, max_elem)
!print*,idx_i,idx_j,max_elem !print*,idx_i,idx_j,max_elem
max_elem = dabs(max_elem) max_elem = dabs(max_elem)
deallocate(cos4theta, sin4theta) deallocate(cos4theta, sin4theta)
deallocate(A,B,beta,gamma) deallocate(A,B,beta,gamma)
end end
subroutine theta_PM(l, n, m_x, max_elem) subroutine theta_PM(l, n, m_x, max_elem)
include 'pi.h' include 'pi.h'
BEGIN_DOC BEGIN_DOC
@ -1507,18 +1507,18 @@ subroutine theta_PM(l, n, m_x, max_elem)
! Loop over the MOs of the a given mo_class to compute <i|P_a|j> ! Loop over the MOs of the a given mo_class to compute <i|P_a|j>
do tmp_j = 1, n do tmp_j = 1, n
j = l(tmp_j) j = l(tmp_j)
do tmp_i = 1, n do tmp_i = 1, n
i = l(tmp_i) i = l(tmp_i)
do rho = 1, ao_num ! loop over all the AOs do rho = 1, ao_num ! loop over all the AOs
do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
mu = nucl_aos(a,b) ! AO centered on atom a mu = nucl_aos(a,b) ! AO centered on atom a
Pa(tmp_i,tmp_j) = Pa(tmp_i,tmp_j) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) & Pa(tmp_i,tmp_j) = Pa(tmp_i,tmp_j) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j)) + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
enddo enddo
enddo enddo
enddo enddo
enddo enddo
@ -1528,7 +1528,7 @@ subroutine theta_PM(l, n, m_x, max_elem)
Aij(i,j) = Aij(i,j) + Pa(i,j)**2 - 0.25d0 * (Pa(i,i) - Pa(j,j))**2 Aij(i,j) = Aij(i,j) + Pa(i,j)**2 - 0.25d0 * (Pa(i,i) - Pa(j,j))**2
enddo enddo
enddo enddo
! B ! B
do j = 1, n do j = 1, n
do i = 1, n do i = 1, n
@ -1581,7 +1581,7 @@ end
! Spatial extent ! Spatial extent
! The spatial extent of an orbital $i$ is computed as ! The spatial extent of an orbital $i$ is computed as
! \begin{align*} ! \begin{align*}
! \sum_{\lambda=x,y,z}\sqrt{<i|\lambda^2|i> - <i|\lambda|i>^2} ! \sum_{\lambda=x,y,z}\sqrt{<i|\lambda^2|i> - <i|\lambda|i>^2}
! \end{align*} ! \end{align*}
@ -1596,14 +1596,14 @@ subroutine compute_spatial_extent(spatial_extent)
BEGIN_DOC BEGIN_DOC
! Compute the spatial extent of the MOs ! Compute the spatial extent of the MOs
END_DOC END_DOC
double precision, intent(out) :: spatial_extent(mo_num) double precision, intent(out) :: spatial_extent(mo_num)
double precision :: average_core, average_act, average_inact, average_virt double precision :: average_core, average_act, average_inact, average_virt
double precision :: std_var_core, std_var_act, std_var_inact, std_var_virt double precision :: std_var_core, std_var_act, std_var_inact, std_var_virt
integer :: i,j,k,l integer :: i,j,k,l
spatial_extent = 0d0 spatial_extent = 0d0
do i = 1, mo_num do i = 1, mo_num
spatial_extent(i) = mo_spread_x(i,i) - mo_dipole_x(i,i)**2 spatial_extent(i) = mo_spread_x(i,i) - mo_dipole_x(i,i)**2
enddo enddo
@ -1631,7 +1631,7 @@ subroutine compute_spatial_extent(spatial_extent)
call compute_average_sp_ext(spatial_extent, list_act, dim_list_act_orb, average_act) call compute_average_sp_ext(spatial_extent, list_act, dim_list_act_orb, average_act)
call compute_std_var_sp_ext(spatial_extent, list_act, dim_list_act_orb, average_act, std_var_act) call compute_std_var_sp_ext(spatial_extent, list_act, dim_list_act_orb, average_act, std_var_act)
endif endif
average_inact = 0d0 average_inact = 0d0
std_var_inact = 0d0 std_var_inact = 0d0
if (dim_list_inact_orb >= 2) then if (dim_list_inact_orb >= 2) then
@ -1661,12 +1661,12 @@ subroutine compute_spatial_extent(spatial_extent)
print*, 'virt:', dim_list_virt_orb print*, 'virt:', dim_list_virt_orb
print*, 'mo_num:', mo_num print*, 'mo_num:', mo_num
print*,'' print*,''
print*,'-- Core MOs --' print*,'-- Core MOs --'
print*,'Average:', average_core print*,'Average:', average_core
print*,'Std var:', std_var_core print*,'Std var:', std_var_core
print*,'' print*,''
print*,'-- Active MOs --' print*,'-- Active MOs --'
print*,'Average:', average_act print*,'Average:', average_act
print*,'Std var:', std_var_act print*,'Std var:', std_var_act
@ -1692,7 +1692,7 @@ end
subroutine compute_average_sp_ext(spatial_extent, list, list_size, average) subroutine compute_average_sp_ext(spatial_extent, list, list_size, average)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute the average spatial extent of the MOs ! Compute the average spatial extent of the MOs
END_DOC END_DOC
@ -1701,7 +1701,7 @@ subroutine compute_average_sp_ext(spatial_extent, list, list_size, average)
double precision, intent(in) :: spatial_extent(mo_num) double precision, intent(in) :: spatial_extent(mo_num)
double precision, intent(out) :: average double precision, intent(out) :: average
integer :: i, tmp_i integer :: i, tmp_i
average = 0d0 average = 0d0
do tmp_i = 1, list_size do tmp_i = 1, list_size
i = list(tmp_i) i = list(tmp_i)
@ -1732,7 +1732,7 @@ subroutine compute_std_var_sp_ext(spatial_extent, list, list_size, average, std_
i = list(tmp_i) i = list(tmp_i)
std_var = std_var + (spatial_extent(i) - average)**2 std_var = std_var + (spatial_extent(i) - average)**2
enddo enddo
std_var = dsqrt(1d0/DBLE(list_size) * std_var) std_var = dsqrt(1d0/DBLE(list_size) * std_var)
end end
@ -1779,7 +1779,7 @@ subroutine apply_pre_rotation()
enddo enddo
enddo enddo
endif endif
! Pre rotation for active MOs ! Pre rotation for active MOs
if (dim_list_act_orb >= 2) then if (dim_list_act_orb >= 2) then
do tmp_j = 1, dim_list_act_orb do tmp_j = 1, dim_list_act_orb
@ -1796,7 +1796,7 @@ subroutine apply_pre_rotation()
enddo enddo
enddo enddo
endif endif
! Pre rotation for inactive MOs ! Pre rotation for inactive MOs
if (dim_list_inact_orb >= 2) then if (dim_list_inact_orb >= 2) then
do tmp_j = 1, dim_list_inact_orb do tmp_j = 1, dim_list_inact_orb
@ -1813,7 +1813,7 @@ subroutine apply_pre_rotation()
enddo enddo
enddo enddo
endif endif
! Pre rotation for virtual MOs ! Pre rotation for virtual MOs
if (dim_list_virt_orb >= 2) then if (dim_list_virt_orb >= 2) then
do tmp_j = 1, dim_list_virt_orb do tmp_j = 1, dim_list_virt_orb
@ -1830,21 +1830,21 @@ subroutine apply_pre_rotation()
enddo enddo
enddo enddo
endif endif
! Nothing for deleted ones ! Nothing for deleted ones
! Compute pre rotation matrix from pre_rot ! Compute pre rotation matrix from pre_rot
call rotation_matrix(pre_rot,mo_num,R,mo_num,mo_num,info,enforce_step_cancellation) call rotation_matrix(pre_rot,mo_num,R,mo_num,mo_num,info,enforce_step_cancellation)
if (enforce_step_cancellation) then if (enforce_step_cancellation) then
print*, 'Cancellation of the pre rotation, too big error in the rotation matrix' print*, 'Cancellation of the pre rotation, too big error in the rotation matrix'
print*, 'Reduce the angle for the pre rotation, abort' print*, 'Reduce the angle for the pre rotation, abort'
call abort call abort
endif endif
! New Mos (we don't car eabout the previous MOs prev_mos) ! New Mos (we don't car eabout the previous MOs prev_mos)
call apply_mo_rotation(R,prev_mos) call apply_mo_rotation(R,prev_mos)
! Update the things related to mo_coef ! Update the things related to mo_coef
TOUCH mo_coef TOUCH mo_coef
call save_mos call save_mos
@ -1885,19 +1885,19 @@ subroutine x_tmp_orb_loc_v2(tmp_n, tmp_list_size, tmp_list, v_grad, H,tmp_x, tmp
! min element in the hessian ! min element in the hessian
if (lambda < 0d0) then if (lambda < 0d0) then
lambda = -lambda + 1d-6 lambda = -lambda + 1d-6
endif endif
print*, 'lambda', lambda print*, 'lambda', lambda
! Good ! Good
do tmp_k = 1, tmp_n do tmp_k = 1, tmp_n
if (ABS(H(tmp_k,tmp_k)) > 1d-6) then if (ABS(H(tmp_k,tmp_k)) > 1d-6) then
tmp_x(tmp_k) = - 1d0/(ABS(H(tmp_k,tmp_k))+lambda) * v_grad(tmp_k)!(-v_grad(tmp_k)) tmp_x(tmp_k) = - 1d0/(ABS(H(tmp_k,tmp_k))+lambda) * v_grad(tmp_k)!(-v_grad(tmp_k))
!x(tmp_k) = - 1d0/(ABS(H(tmp_k,tmp_k))+lambda) * (-v_grad(tmp_k)) !x(tmp_k) = - 1d0/(ABS(H(tmp_k,tmp_k))+lambda) * (-v_grad(tmp_k))
endif endif
enddo enddo
! 1D tmp -> 2D tmp ! 1D tmp -> 2D tmp
tmp_m_x = 0d0 tmp_m_x = 0d0
do tmp_j = 1, tmp_list_size - 1 do tmp_j = 1, tmp_list_size - 1
do tmp_i = tmp_j + 1, tmp_list_size do tmp_i = tmp_j + 1, tmp_list_size
@ -1909,7 +1909,7 @@ subroutine x_tmp_orb_loc_v2(tmp_n, tmp_list_size, tmp_list, v_grad, H,tmp_x, tmp
! Antisym ! Antisym
do tmp_i = 1, tmp_list_size - 1 do tmp_i = 1, tmp_list_size - 1
do tmp_j = tmp_i + 1, tmp_list_size do tmp_j = tmp_i + 1, tmp_list_size
tmp_m_x(tmp_i,tmp_j) = - tmp_m_x(tmp_j,tmp_i) tmp_m_x(tmp_i,tmp_j) = - tmp_m_x(tmp_j,tmp_i)
enddo enddo
enddo enddo
@ -1947,16 +1947,17 @@ subroutine ao_to_mo_no_sym(A_ao,LDA_ao,A_mo,LDA_mo)
end end
subroutine run_sort_by_fock_energies() subroutine run_sort_by_fock_energies()
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Saves the current MOs ordered by diagonal element of the Fock operator. ! Saves the current MOs ordered by diagonal element of the Fock operator.
END_DOC END_DOC
integer :: i,j,k,l,tmp_i,tmp_k,tmp_list_size integer :: i,j,k,l,tmp_i,tmp_k,tmp_list_size
integer, allocatable :: iorder(:) integer, allocatable :: iorder(:)
double precision, allocatable :: fock_energies_tmp(:), tmp_mo_coef(:,:), tmp_list(:) double precision, allocatable :: fock_energies_tmp(:), tmp_mo_coef(:,:)
integer, allocatable :: tmp_list(:)
! allocate(iorder(mo_num), fock_energies_tmp(mo_num), new_mo_coef(ao_num, mo_num)) ! allocate(iorder(mo_num), fock_energies_tmp(mo_num), new_mo_coef(ao_num, mo_num))
! !
@ -2006,7 +2007,7 @@ subroutine run_sort_by_fock_energies()
else else
tmp_list = list_virt tmp_list = list_virt
endif endif
print*,'MO class: ',trim(mo_class(tmp_list(1))) print*,'MO class: ', trim(mo_class(tmp_list(1)))
allocate(iorder(tmp_list_size), fock_energies_tmp(tmp_list_size), tmp_mo_coef(ao_num,tmp_list_size)) allocate(iorder(tmp_list_size), fock_energies_tmp(tmp_list_size), tmp_mo_coef(ao_num,tmp_list_size))
!print*,'MOs before sorting them by f_p^p energies:' !print*,'MOs before sorting them by f_p^p energies:'
@ -2016,9 +2017,9 @@ subroutine run_sort_by_fock_energies()
iorder(i) = i iorder(i) = i
!print*, tmp_i, fock_energies_tmp(i) !print*, tmp_i, fock_energies_tmp(i)
enddo enddo
call dsort(fock_energies_tmp, iorder, tmp_list_size) call dsort(fock_energies_tmp, iorder, tmp_list_size)
print*,'MOs after sorting them by f_p^p energies:' print*,'MOs after sorting them by f_p^p energies:'
do i = 1, tmp_list_size do i = 1, tmp_list_size
k = iorder(i) k = iorder(i)
@ -2042,7 +2043,7 @@ subroutine run_sort_by_fock_energies()
print*,'HF energy:', HF_energy print*,'HF energy:', HF_energy
endif endif
print*,'' print*,''
deallocate(iorder, fock_energies_tmp, tmp_list, tmp_mo_coef) deallocate(iorder, fock_energies_tmp, tmp_list, tmp_mo_coef)
endif endif
@ -2050,5 +2051,5 @@ subroutine run_sort_by_fock_energies()
touch mo_coef touch mo_coef
call save_mos call save_mos
end end