* Orbital localization Molecular orbitals localization ** Doc The program localizes the orbitals in function of their mo_class: - core MOs - inactive MOs - active MOs - virtual MOs - deleted MOs -> no orbital localization Core MOs are localized with core MOs, inactives MOs are localized with inactives MOs and so on. But deleted orbitals are not localized. WARNING: - The user MUST SPECIFY THE MO CLASSES, otherwise if default mo class is false the localization will be done for all the orbitals between them, so the occupied and virtual MOs will be combined together which is clearly not what we want to do. If default lpmo class is true the localization will be done for the core, occupied and virtual orbitals, but pay attention the mo_class are not deleted after... - The mo class is not important (except "deleted") because it is not link to the kind of MOs for CASSCF or CIPSI. It is just a way to separate the MOs in order to localize them separetely, for example to separate the core MOs, the occupied MOs and the virtuals MOs. - The user MUST CHANGE THE MO CLASSES AFTER THE LOCALIZATION in order to have the right mo class for his next calculation... For more information on the mo_class: lpqp set_mo_class -h *** Foster-Boys localization Foster-Boys localization: - cite Foster Boys, S. F., 1960, Rev. Mod. Phys. 32, 296. DOI:https://doi.org/10.1103/RevModPhys.32.300 Boys, S. F., 1966, in Quantum Theory of Atoms, Molecules, and the Solid State, edited by P.-O. Löwdin (Academic Press, New York), p. 253. Daniel A. Kleier, Thomas A. Halgren, John H. Hall Jr., and William N. Lipscomb, J. Chem. Phys. 61, 3905 (1974) doi: 10.1063/1.1681683 Høyvik, I.-M., Jansik, B., Jørgensen, P., J. Comput. Chem. 2013, 34, 1456– 1462. DOI: 10.1002/jcc.23281 Høyvik, I.-M., Jansik, B., Jørgensen, P., J. Chem. Theory Comput. 2012, 8, 9, 3137–3146 DOI: https://doi.org/10.1021/ct300473g Høyvik, I.-M., Jansik, B., Jørgensen, P., J. Chem. Phys. 137, 224114 (2012) DOI: https://doi.org/10.1063/1.4769866 Nicola Marzari, Arash A. Mostofi, Jonathan R. Yates, Ivo Souza, and David Vanderbilt Rev. Mod. Phys. 84, 1419 https://doi.org/10.1103/RevModPhys.84.1419 The Foster-Boys localization is a method to generate localized MOs (LMOs) by minimizing the Foster-Boys criterion: $$ C_{FB} = \sum_{i=1}^N \left[ < \phi_i | r^2 | \phi_i > - < \phi_i | r | \phi_i >^2 \right] $$. In fact it is equivalent to maximise $$ C_2 = \sum_{i>j, \ i=1}^N \left[ < \phi_i | r | \phi_i > - < \phi_j | r | \phi_j > \left]^2$$ or $$ C_3 = \sum_{i=1}^N \left[ < \phi_i | r | \phi_i > \right]^2.$$ Noting $$A_{ii} = < \phi_i | r^2 | \phi_i > $$ $$B_{ii} = < \phi_i | r | \phi_i > $$ $$ \beta = (B_{pp} - B_{qq})^2 - 4 B_{pq}^2 $$ $$ \gamma = 4 B_{pq} (B_{pp} - B_{qq}) $$ \begin{align*} C_{FB}(\theta) &= \sum_{i=1}^N \left[ A_{ii} - B_{ii}^2 \right] \\ &- \left[ A_{pp} - B_{pp}^2 + A_{qq} - B_{qq}^2 \right] \\ &+ \left[ A_{pp} + A_{qq} - B_{pp}^2 - B_{qq}^2 + \frac{1}{4} [(1-\cos(4\theta) \beta + \sin(4\theta) \gamma] \right] \\ &= C_1(\theta=0) + \frac{1}{4} [(1-\cos(4\theta)) \beta + \sin(4\theta) \gamma] \end{align*} The derivatives are: \begin{align*} \frac{\partial C_{FB}(\theta)}{\partial \theta} = \beta \sin(4\theta) + \gamma \cos(4 \theta) \end{align*} \begin{align*} \frac{\partial^2 C_{FB}(\theta)}{\partial \theta^2} = 4 \beta \cos(4\theta) - 4 \gamma \sin(4 \theta) \end{align*} Similarly: \begin{align*} C_3(\theta) &= \sum_{i=1}^N [B_{ii}^2] \\ &- B_{pp}^2 - B_{qq}^2 \\ &+ B_{pp}^2 + B_{qq}^2 - \frac{1}{4} [(1-\cos(4\theta) \beta + \sin(4\theta) \gamma] \\ &= C_3(\theta=0) - \frac{1}{4} [(1-\cos(4\theta)) \beta + \sin(4\theta) \gamma] \end{align*} The derivatives are: \begin{align*} \frac{\partial C_3(\theta)}{\partial \theta} = - \beta \sin(4\theta) - \gamma \cos(4 \theta) \end{align*} \begin{align*} \frac{\partial^2 C_3(\theta)}{\partial \theta^2} = - 4 \beta \cos(4\theta) + 4 \gamma \sin(4 \theta) \end{align*} And since we compute the derivatives around $\theta = 0$ (around the actual position) we have: \begin{align*} \left. \frac{\partial{C_{FB}(\theta)}}{\partial \theta}\right|_{\theta=0} = \gamma \end{align*} \begin{align*} \left. \frac{\partial^2 C_{FB}(\theta)}{\partial \theta^2}\right|_{\theta=0} = 4 \beta \end{align*} Locality of the orbitals: - cite Hoyvik As the Foster-Boys method tries to minimize the sum of the second moment MO spread, the locality of each MO can be expressed as the second moment of the MO spread. For the MO i, the locality criterion is \begin{align*} \sigma_i &= \sqrt{ - ^2} \\ &= \sqrt{ - ^2 + - ^2 + - ^2} \end{align*} *** Pipek-Mezey localization -cite pipek mezey 1989 J. Pipek, P. G. Mezey, J. Chem. Phys. 90, 4916 (1989) DOI: 10.1063/1.456588 Foster-Boys localization does not preserve the $\sigma - \pi$ separation of the MOs, it leads to "banana" orbitals. The Pipek-Mezey localization normally preserves this separation. The optimum functional $\mathcal{P}$ is obtained for the maximum of $D^{-1}$ \begin{align*} \mathcal{P} = \sum_{i=1}^n \sum_{A=1}^N \left[ \right]^2 \end{align*} As for the Foster Boys localization, the change in the functional for the rotation of two MOs can be obtained using very similar terms \begin{align*} \beta_{st}^{PM} = \sum_{A=1}^N \left( ^2 - \frac{1}{4} \left[ - \right]^2 \right) \end{align*} \begin{align*} \gamma_{st}^{PM} = \sum_{A=1}^N \left[ - \right] \end{align*} The matrix element of the operator $P_A$ are obtained using \begin{align*} <\rho | \tilde{\mu}> = \delta_{\rho \mu} \end{align*} which leads to \begin{align*} = \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*} $\sum_{\rho}$ -> sum over all the AOs $\sum_{\mu \in A}$ -> sum over the AOs which belongs to atom A $c^t$ -> expansion coefficient of orbital |t> So similarly the first and second derivatives are \begin{align*} \left. \frac{\partial \mathcal{P} (\theta)}{\partial \theta} \right|_{\theta=0}= \gamma^{PM} \end{align*} \begin{align*} \left. \frac{\partial^2 \mathcal{P} (\theta)}{\partial \theta^2}\right|_{\theta=0} = 4 \beta^{PM} \end{align*} ** Localization procedure Localization procedure: To do the localization we compute the gradient and the diagonal hessian of the Foster-Boys criterion with respect to the MO rotations and we minimize it with the Newton method. In order to avoid the problem of starting on a saddle point, the localization procedure starts by giving a little kick in the MOs, by putting "kick in mos" true, in order to break the symmetry and escape from a possible saddle point. In order to speed up the iteration we compute the gradient, the diagonal hessian and the step in temporary matrices of the size (number MOs in mo class by number MOs in mo class) ** Remarks Variables: The indexes i and j refere to the positions of the elements in the "full space", i.e., the arrays containing elements for all the MOs, but the indexes tmp_i and tmp_j to the positions of the elements in the "reduced space/subspace", i.e., the arrays containing elements for a restricted number of MOs. Example: The gradient for the localization of the core MOs can be expressed as a vector of length mo_num*(mo_num-1)/2 with only n_core_orb*(n_core_orb-1)/2 non zero elements, so it is more relevant to use a vector of size n_act_orb*(n_core_orb-1)/2. So here the gradient is a vector of size tmp_list_size*(tmp_list_size)/2 where tmp_list_size is the number of MOs is the corresponding mo class. The same thing happened for the hessian, the matrix containing the step and the rotation matrix, which are tmp_list_size by tmp_list_size matrices. Ex gradient for 4 core orbitales: \begin{align*} \begin{pmatrix} 0 & -a & -b & -d & \hdots & 0 \\ a & 0 & -c & -e & \hdots & 0 \\ b & c & 0 & -f & \hdots & 0 \\ d & e & f & 0 & \hdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \hdots & 0 \\ \end{pmatrix} \Rightarrow \begin{pmatrix} a \\ b \\ c \\ e \\ f \\ 0 \\ \vdots \\ 0 \\ \end{pmatrix} \end{align*} \begin{align*} \begin{pmatrix} 0 & -a & -b & -d & \hdots & 0 \\ a & 0 & -c & -e & \hdots & 0 \\ b & c & 0 & -f & \hdots & 0 \\ d & e & f & 0 & \hdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \hdots & 0 \\ \end{pmatrix} \Rightarrow \begin{pmatrix} 0 & -a & -b & -d \\ a & 0 & -c & -e \\ b & c & 0 & -f \\ d & e & f & 0 \\ \end{pmatrix} \Rightarrow \begin{pmatrix} a \\ b \\ c \\ e \\ f \\ \end{pmatrix} \end{align*} The same thing can be done if indexes of the orbitales are not consecutives since it's done with lists of MOs: \begin{align*} \begin{pmatrix} 0 & -a & 0 & -b & -d & \hdots & 0 \\ a & 0 & 0 & -c & -e & \hdots & 0 \\ 0 & 0 & 0 & 0 & 0 & \hdots & 0 \\ b & c & 0 & 0 & -f & \hdots & 0 \\ d & e & 0 & f & 0 & \hdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \hdots & 0 \\ \end{pmatrix} \Rightarrow \begin{pmatrix} 0 & -a & -b & -d \\ a & 0 & -c & -e \\ b & c & 0 & -f \\ d & e & f & 0 \\ \end{pmatrix} \Rightarrow \begin{pmatrix} a \\ b \\ c \\ e \\ f \\ \end{pmatrix} \end{align*} The dipoles are updated using the "ao to mo" subroutine without the "restore symmetry" which is actually in N^4 but can be rewrite in N^2 log(N^2). The bottleneck of the program is normally N^3 with the matrix multiplications/diagonalizations. The use of the full hessian can be an improvement but it will scale in N^4... ** Program #+BEGIN_SRC f90 org :tangle localization.irp.f program localization implicit none call set_classes_loc call run_localization call unset_classes_loc end #+END_SRC Variables: | pre_rot(mo_num, mo_num) | double precision | Matrix for the pre rotation | | R(mo_num,mo_num) | double precision | Rotation matrix | | tmp_R(:,:) | double precision | Rottation matrix in a subsapce | | prev_mos(ao_num, mo_num) | double precision | Previous mo_coef | | spatial_extent(mo_num) | double precision | Spatial extent of the orbitals | | criterion | double precision | Localization criterion | | prev_criterion | double precision | Previous criterion | | criterion_model | double precision | Estimated next criterion | | rho | double precision | Ratio to measure the agreement between the model | | | | and the reality | | delta | double precision | Radisu of the trust region | | norm_grad | double precision | Norm of the gradient | | info | integer | for dsyev from Lapack | | max_elem | double precision | maximal element in the gradient | | v_grad(:) | double precision | Gradient | | H(:,:) | double precision | Hessian (diagonal) | | e_val(:) | double precision | Eigenvalues of the hessian | | W(:,:) | double precision | Eigenvectors of the hessian | | tmp_x(:) | double precision | Step in 1D (in a subaspace) | | tmp_m_x(:,:) | double precision | Step in 2D (in a subaspace) | | tmp_list(:) | double precision | List of MOs in a mo_class | | i,j,k | integer | Indexes in the full MO space | | tmp_i, tmp_j, tmp_k | integer | Indexes in a subspace | | l | integer | Index for the mo_class | | key(:) | integer | Key to sort the eigenvalues of the hessian | | nb_iter | integer | Number of iterations | | must_exit | logical | To exit the trust region loop | | cancel_step | logical | To cancel a step | | not_*converged | logical | To localize the different mo classes | | t* | double precision | To measure the time | | n | integer | mo_num*(mo_num-1)/2, number of orbital parameters | | tmp_n | integer | dim_subspace*(dim_subspace-1)/2 | | | | Number of dimension in the subspace | Variables in qp_edit for the localization: | localization_method | | localization_max_nb_iter | | default_mo_class | | thresh_loc_max_elem_grad | | kick_in_mos | | angle_pre_rot | + all the variables for the trust region Cf. qp_edit orbital optimization #+BEGIN_SRC f90 :comments org :tangle localization.irp.f subroutine run_localization include 'pi.h' BEGIN_DOC ! Orbital localization END_DOC implicit none ! Variables double precision, allocatable :: pre_rot(:,:), R(:,:) double precision, allocatable :: prev_mos(:,:), spatial_extent(:), tmp_R(:,:) double precision :: criterion, norm_grad integer :: i,j,k,l,p, tmp_i, tmp_j, tmp_k integer :: info integer :: n, tmp_n, tmp_list_size double precision, allocatable :: v_grad(:), H(:), tmp_m_x(:,:), tmp_x(:),W(:),e_val(:) double precision :: max_elem, t1, t2, t3, t4, t5, t6 integer, allocatable :: tmp_list(:), key(:) double precision :: prev_criterion, rho, delta, criterion_model integer :: nb_iter, nb_sub_iter logical :: not_converged, not_core_converged logical :: not_act_converged, not_inact_converged, not_virt_converged logical :: use_trust_region, must_exit, cancel_step,enforce_step_cancellation n = mo_num*(mo_num-1)/2 ! Allocation allocate(spatial_extent(mo_num)) allocate(pre_rot(mo_num, mo_num), R(mo_num, mo_num)) allocate(prev_mos(ao_num, mo_num)) ! Locality before the localization call compute_spatial_extent(spatial_extent) ! Choice of the method print*,'' print*,'Localization method:',localization_method if (localization_method == 'boys') then print*,'Foster-Boys localization' elseif (localization_method == 'pipek') then print*,'Pipek-Mezey localization' else print*,'Unknown localization_method, please select boys or pipek' call abort endif print*,'' ! Localization criterion (FB, PM, ...) for each mo_class print*,'### Before the pre rotation' ! Debug if (debug_hf) then print*,'HF energy:', HF_energy endif do l = 1, 4 if (l==1) then ! core tmp_list_size = dim_list_core_orb elseif (l==2) then ! act tmp_list_size = dim_list_act_orb elseif (l==3) then ! inact tmp_list_size = dim_list_inact_orb else ! virt tmp_list_size = dim_list_virt_orb endif ! Allocation tmp array allocate(tmp_list(tmp_list_size)) ! To give the list of MOs in a mo_class if (l==1) then ! core tmp_list = list_core elseif (l==2) then tmp_list = list_act elseif (l==3) then tmp_list = list_inact else tmp_list = list_virt endif if (tmp_list_size >= 2) then call criterion_localization(tmp_list_size, tmp_list,criterion) print*,'Criterion:', criterion, mo_class(tmp_list(1)) endif deallocate(tmp_list) enddo ! Debug !print*,'HF', HF_energy #+END_SRC ** Loc #+BEGIN_SRC f90 :comments org :tangle localization.irp.f ! Pre rotation, to give a little kick in the MOs call apply_pre_rotation() ! Criterion after the pre rotation ! Localization criterion (FB, PM, ...) for each mo_class print*,'### After the pre rotation' ! Debug if (debug_hf) then touch mo_coef print*,'HF energy:', HF_energy endif do l = 1, 4 if (l==1) then ! core tmp_list_size = dim_list_core_orb elseif (l==2) then ! act tmp_list_size = dim_list_act_orb elseif (l==3) then ! inact tmp_list_size = dim_list_inact_orb else ! virt tmp_list_size = dim_list_virt_orb endif if (tmp_list_size >= 2) then ! Allocation tmp array allocate(tmp_list(tmp_list_size)) ! To give the list of MOs in a mo_class if (l==1) then ! core tmp_list = list_core elseif (l==2) then tmp_list = list_act elseif (l==3) then tmp_list = list_inact else tmp_list = list_virt endif call criterion_localization(tmp_list_size, tmp_list,criterion) print*,'Criterion:', criterion, trim(mo_class(tmp_list(1))) deallocate(tmp_list) endif enddo ! Debug !print*,'HF', HF_energy print*,'' print*,'========================' print*,' Orbital localization' print*,'========================' print*,'' !Initialization not_converged = .TRUE. ! To do the localization only if there is at least 2 MOs if (dim_list_core_orb >= 2) then not_core_converged = .TRUE. else not_core_converged = .FALSE. endif if (dim_list_act_orb >= 2) then not_act_converged = .TRUE. else not_act_converged = .FALSE. endif if (dim_list_inact_orb >= 2) then not_inact_converged = .TRUE. else not_inact_converged = .FALSE. endif if (dim_list_virt_orb >= 2) then not_virt_converged = .TRUE. else not_virt_converged = .FALSE. endif ! Loop over the mo_classes do l = 1, 4 if (l==1) then ! core not_converged = not_core_converged tmp_list_size = dim_list_core_orb elseif (l==2) then ! act not_converged = not_act_converged tmp_list_size = dim_list_act_orb elseif (l==3) then ! inact not_converged = not_inact_converged tmp_list_size = dim_list_inact_orb else ! virt not_converged = not_virt_converged tmp_list_size = dim_list_virt_orb endif ! Next iteration if converged = true if (.not. not_converged) then cycle endif ! Allocation tmp array allocate(tmp_list(tmp_list_size)) ! To give the list of MOs in a mo_class if (l==1) then ! core tmp_list = list_core elseif (l==2) then tmp_list = list_act elseif (l==3) then tmp_list = list_inact else tmp_list = list_virt endif ! Display if (not_converged) then print*,'' print*,'###', trim(mo_class(tmp_list(1))), 'MOs ###' print*,'' endif ! Size for the 2D -> 1D transformation tmp_n = tmp_list_size * (tmp_list_size - 1)/2 ! Without hessian + trust region if (.not. localization_use_hessian) then ! Allocation of temporary arrays allocate(v_grad(tmp_n), tmp_m_x(tmp_list_size, tmp_list_size)) allocate(tmp_R(tmp_list_size, tmp_list_size), tmp_x(tmp_n)) ! Criterion call criterion_localization(tmp_list_size, tmp_list, prev_criterion) ! Init nb_iter = 0 delta = 1d0 !Loop do while (not_converged) print*,'' print*,'***********************' print*,'Iteration', nb_iter print*,'***********************' print*,'' ! Angles of rotation call theta_localization(tmp_list, tmp_list_size, tmp_m_x, max_elem) tmp_m_x = - tmp_m_x * delta ! Rotation submatrix call rotation_matrix(tmp_m_x, tmp_list_size, tmp_R, tmp_list_size, tmp_list_size, & info, enforce_step_cancellation) ! To ensure that the rotation matrix is unitary if (enforce_step_cancellation) then print*, 'Step cancellation, too large error in the rotation matrix' delta = delta * 0.5d0 cycle else delta = min(delta * 2d0, 1d0) endif ! Full rotation matrix and application of the rotation call sub_to_full_rotation_matrix(tmp_list_size, tmp_list, tmp_R, R) call apply_mo_rotation(R, prev_mos) ! Update the needed data call update_data_localization() ! New criterion call criterion_localization(tmp_list_size, tmp_list, criterion) print*,'Criterion:', trim(mo_class(tmp_list(1))), nb_iter, criterion print*,'Max elem :', max_elem print*,'Delta :', delta nb_iter = nb_iter + 1 ! Exit if (nb_iter >= localization_max_nb_iter .or. dabs(max_elem) < thresh_loc_max_elem_grad) then not_converged = .False. endif enddo ! Save the changes call update_data_localization() call save_mos() TOUCH mo_coef ! Deallocate deallocate(v_grad, tmp_m_x, tmp_list) deallocate(tmp_R, tmp_x) ! Trust region else ! Allocation of temporary arrays allocate(v_grad(tmp_n), H(tmp_n), tmp_m_x(tmp_list_size, tmp_list_size)) allocate(tmp_R(tmp_list_size, tmp_list_size)) allocate(tmp_x(tmp_n), W(tmp_n), e_val(tmp_n), key(tmp_n)) ! ### Initialization ### delta = 0d0 ! can be deleted (normally) nb_iter = 0 ! Must start at 0 !!! rho = 0.5d0 ! Must be 0.5 ! Compute the criterion before the loop call criterion_localization(tmp_list_size, tmp_list, prev_criterion) ! Loop until the convergence do while (not_converged) print*,'' print*,'***********************' print*,'Iteration', nb_iter print*,'***********************' print*,'' ! Gradient call gradient_localization(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) ! Diagonal hessian call hessian_localization(tmp_n, tmp_list_size, tmp_list, H) ! Diagonalization of the diagonal hessian by hands !call diagonalization_hessian(tmp_n,H,e_val,w) do i = 1, tmp_n e_val(i) = H(i) enddo ! Key list for dsort do i = 1, tmp_n key(i) = i enddo ! Sort of the eigenvalues call dsort(e_val, key, tmp_n) ! Eigenvectors W = 0d0 do i = 1, tmp_n W(i) = dble(key(i)) enddo ! To enter in the loop just after cancel_step = .True. nb_sub_iter = 0 ! Loop to reduce the trust radius until the criterion decreases and rho >= thresh_rho do while (cancel_step) print*,'-----------------------------' print*, mo_class(tmp_list(1)) print*,'Iteration:', nb_iter print*,'Sub iteration:', nb_sub_iter print*,'Max elem grad:', max_elem print*,'-----------------------------' ! Hessian,gradient,Criterion -> x call trust_region_step_w_expected_e(tmp_n,1, H, W, e_val, v_grad, prev_criterion, & rho, nb_iter, delta, criterion_model, tmp_x, must_exit) ! Internal loop exit condition if (must_exit) then print*,'trust_region_step_w_expected_e sent: Exit' exit endif ! 1D tmp -> 2D tmp call vec_to_mat_v2(tmp_n, tmp_list_size, tmp_x, tmp_m_x) ! Rotation submatrix (square matrix tmp_list_size by tmp_list_size) call rotation_matrix(tmp_m_x, tmp_list_size, tmp_R, tmp_list_size, tmp_list_size, & info, enforce_step_cancellation) if (enforce_step_cancellation) then print*, 'Step cancellation, too large error in the rotation matrix' rho = 0d0 cycle endif ! tmp_R to R, subspace to full space call sub_to_full_rotation_matrix(tmp_list_size, tmp_list, tmp_R, R) ! Rotation of the MOs call apply_mo_rotation(R, prev_mos) ! Update the things related to mo_coef call update_data_localization() ! Update the criterion call criterion_localization(tmp_list_size, tmp_list, criterion) print*,'Criterion:', trim(mo_class(tmp_list(1))), nb_iter, criterion ! Criterion -> step accepted or rejected call trust_region_is_step_cancelled(nb_iter, prev_criterion, criterion, & criterion_model, rho, cancel_step) ! Cancellation of the step, previous MOs if (cancel_step) then mo_coef = prev_mos endif nb_sub_iter = nb_sub_iter + 1 enddo !call save_mos() !### depend of the time for 1 iteration ! To exit the external loop if must_exti = .True. if (must_exit) then exit endif ! Step accepted, nb iteration + 1 nb_iter = nb_iter + 1 ! External loop exit conditions if (DABS(max_elem) < thresh_loc_max_elem_grad) then not_converged = .False. endif if (nb_iter > localization_max_nb_iter) then not_converged = .False. endif enddo ! Deallocation of temporary arrays deallocate(v_grad, H, tmp_m_x, tmp_R, tmp_list, tmp_x, W, e_val, key) ! Save the MOs call save_mos() TOUCH mo_coef ! Debug if (debug_hf) then touch mo_coef print*,'HF energy:', HF_energy endif endif enddo ! Seems unecessary TOUCH mo_coef ! To sort the MOs using the diagonal elements of the Fock matrix if (sort_mos_by_e) then call run_sort_by_fock_energies() endif ! Debug if (debug_hf) then touch mo_coef print*,'HF energy:', HF_energy endif ! Locality after the localization call compute_spatial_extent(spatial_extent) end #+END_SRC ** Gathering Gradient/hessian/criterion for the localization: They are chosen in function of the localization method Gradient: qp_edit : | localization_method | method for the localization | Input: | 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(tmp_list_size) | integer | MOs in the mo_class | Output: | v_grad(tmp_n) | double precision | Gradient in the subspace | | max_elem | double precision | Maximal element in the gradient | | norm_grad | double precision | Norm of the gradient | #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine gradient_localization(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) include 'pi.h' implicit none BEGIN_DOC ! Compute the gradient of the chosen localization method END_DOC 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 if (localization_method == 'boys') then call gradient_FB_omp(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) !call gradient_FB(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) elseif (localization_method== 'pipek') then call gradient_PM(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) else print*,'Unkown method:'//localization_method call abort endif end #+END_SRC Hessian: Output: | H(tmp_n,tmp_n) | double precision | Gradient in the subspace | | max_elem | double precision | Maximal element in the gradient | | norm_grad | double precision | Norm of the gradient | #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine hessian_localization(tmp_n, tmp_list_size, tmp_list, H) include 'pi.h' implicit none BEGIN_DOC ! Compute the diagonal hessian of the chosen localization method END_DOC integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) double precision, intent(out) :: H(tmp_n) if (localization_method == 'boys') then call hessian_FB_omp(tmp_n, tmp_list_size, tmp_list, H) !call hessian_FB(tmp_n, tmp_list_size, tmp_list, H) ! non OMP for debugging elseif (localization_method == 'pipek') then call hessian_PM(tmp_n, tmp_list_size, tmp_list, H) else print*,'Unkown method: '//localization_method call abort endif end #+END_SRC Criterion: Output: | criterion | double precision | Criterion for the orbital localization | #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine criterion_localization(tmp_list_size, tmp_list,criterion) include 'pi.h' implicit none BEGIN_DOC ! Compute the localization criterion of the chosen localization method END_DOC integer, intent(in) :: tmp_list_size, tmp_list(tmp_list_size) double precision, intent(out) :: criterion if (localization_method == 'boys') then call criterion_FB(tmp_list_size, tmp_list, criterion) elseif (localization_method == 'pipek') then !call criterion_PM(tmp_list_size, tmp_list,criterion) call criterion_PM_v3(tmp_list_size, tmp_list, criterion) else print*,'Unkown method: '//localization_method call abort endif end #+END_SRC Subroutine to update the datas needed for the localization #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine update_data_localization() include 'pi.h' implicit none if (localization_method == 'boys') then ! Update the dipoles call ao_to_mo_no_sym(ao_dipole_x, ao_num, mo_dipole_x, mo_num) call ao_to_mo_no_sym(ao_dipole_y, ao_num, mo_dipole_y, mo_num) call ao_to_mo_no_sym(ao_dipole_z, ao_num, mo_dipole_z, mo_num) elseif (localization_method == 'pipek') then ! Nothing required else print*,'Unkown method: '//localization_method call abort endif end #+END_SRC Angles: Output: | tmp_m_x(tmp_list_size, tmp_list_size) | double precision | Angles for the rotations in the subspace | | max_elem | double precision | Maximal angle | #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine theta_localization(tmp_list, tmp_list_size, tmp_m_x, max_elem) include 'pi.h' implicit none BEGIN_DOC ! Compute the rotation angles between the MOs for the chosen localization method END_DOC 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 if (localization_method == 'boys') then call theta_FB(tmp_list, tmp_list_size, tmp_m_x, max_elem) elseif (localization_method== 'pipek') then call theta_PM(tmp_list, tmp_list_size, tmp_m_x, max_elem) else print*,'Unkown method: '//localization_method call abort endif end #+END_SRC ** Foster-Boys *** Gradient Input: | 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(tmp_list_size) | integer | MOs in the mo_class | Output: | v_grad(tmp_n) | double precision | Gradient in the subspace | | max_elem | double precision | Maximal element in the gradient | | norm_grad | double precision | Norm of the gradient | Internal: | m_grad(tmp_n,tmp_n) | double precision | Gradient in the matrix form | | i,j,k | integer | indexes in the full space | | tmp_i,tmp_j,tmp_k | integer | indexes in the subspace | | t* | double precision | to compute the time | #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine gradient_FB(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) implicit none BEGIN_DOC ! Compute the gradient for the Foster-Boys localization END_DOC 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, allocatable :: m_grad(:,:) integer :: i,j,k,tmp_i,tmp_j,tmp_k double precision :: t1, t2, t3 print*,'' print*,'---gradient_FB---' call wall_time(t1) ! Allocation allocate(m_grad(tmp_list_size, tmp_list_size)) ! Calculation do tmp_j = 1, tmp_list_size j = tmp_list(tmp_j) do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) m_grad(tmp_i,tmp_j) = 4d0 * mo_dipole_x(i,j) * (mo_dipole_x(i,i) - mo_dipole_x(j,j)) & +4d0 * mo_dipole_y(i,j) * (mo_dipole_y(i,i) - mo_dipole_y(j,j)) & +4d0 * mo_dipole_z(i,j) * (mo_dipole_z(i,i) - mo_dipole_z(j,j)) enddo enddo ! 2D -> 1D do tmp_k = 1, tmp_n call vec_to_mat_index(tmp_k,tmp_i,tmp_j) v_grad(tmp_k) = m_grad(tmp_i,tmp_j) enddo ! Maximum element in the gradient max_elem = 0d0 do tmp_k = 1, tmp_n if (ABS(v_grad(tmp_k)) > max_elem) then max_elem = ABS(v_grad(tmp_k)) endif enddo ! Norm of the gradient norm_grad = 0d0 do tmp_k = 1, tmp_n norm_grad = norm_grad + v_grad(tmp_k)**2 enddo norm_grad = dsqrt(norm_grad) print*, 'Maximal element in the gradient:', max_elem print*, 'Norm of the gradient:', norm_grad ! Deallocation deallocate(m_grad) call wall_time(t2) t3 = t2 - t1 print*,'Time in gradient_FB:', t3 print*,'---End gradient_FB---' end subroutine #+END_SRC *** Gradient (OMP) #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine gradient_FB_omp(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) use omp_lib implicit none BEGIN_DOC ! Compute the gradient for the Foster-Boys localization END_DOC 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, allocatable :: m_grad(:,:) integer :: i,j,k,tmp_i,tmp_j,tmp_k double precision :: t1, t2, t3 print*,'' print*,'---gradient_FB_omp---' call wall_time(t1) ! Allocation allocate(m_grad(tmp_list_size, tmp_list_size)) ! Initialization omp call omp_set_max_active_levels(1) !$OMP PARALLEL & !$OMP PRIVATE(i,j,tmp_i,tmp_j,tmp_k) & !$OMP SHARED(tmp_n,tmp_list_size,m_grad,v_grad,mo_dipole_x,mo_dipole_y,mo_dipole_z,tmp_list) & !$OMP DEFAULT(NONE) ! Calculation !$OMP DO do tmp_j = 1, tmp_list_size j = tmp_list(tmp_j) do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) m_grad(tmp_i,tmp_j) = 4d0 * mo_dipole_x(i,j) * (mo_dipole_x(i,i) - mo_dipole_x(j,j)) & +4d0 * mo_dipole_y(i,j) * (mo_dipole_y(i,i) - mo_dipole_y(j,j)) & +4d0 * mo_dipole_z(i,j) * (mo_dipole_z(i,i) - mo_dipole_z(j,j)) enddo enddo !$OMP END DO ! 2D -> 1D !$OMP DO do tmp_k = 1, tmp_n call vec_to_mat_index(tmp_k,tmp_i,tmp_j) v_grad(tmp_k) = m_grad(tmp_i,tmp_j) enddo !$OMP END DO !$OMP END PARALLEL call omp_set_max_active_levels(4) ! Maximum element in the gradient max_elem = 0d0 do tmp_k = 1, tmp_n if (ABS(v_grad(tmp_k)) > max_elem) then max_elem = ABS(v_grad(tmp_k)) endif enddo ! Norm of the gradient norm_grad = 0d0 do tmp_k = 1, tmp_n norm_grad = norm_grad + v_grad(tmp_k)**2 enddo norm_grad = dsqrt(norm_grad) print*, 'Maximal element in the gradient:', max_elem print*, 'Norm of the gradient:', norm_grad ! Deallocation deallocate(m_grad) call wall_time(t2) t3 = t2 - t1 print*,'Time in gradient_FB_omp:', t3 print*,'---End gradient_FB_omp---' end subroutine #+END_SRC *** Hessian Output: | H(tmp_n,tmp_n) | double precision | Gradient in the subspace | | max_elem | double precision | Maximal element in the gradient | | norm_grad | double precision | Norm of the gradient | Internal: Internal: | beta(tmp_n,tmp_n) | double precision | beta in the documentation below to compute the hesian | | i,j,k | integer | indexes in the full space | | tmp_i,tmp_j,tmp_k | integer | indexes in the subspace | | t* | double precision | to compute the time | #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine hessian_FB(tmp_n, tmp_list_size, tmp_list, H) implicit none BEGIN_DOC ! Compute the diagonal hessian for the Foster-Boys localization END_DOC integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) double precision, intent(out) :: H(tmp_n) double precision, allocatable :: beta(:,:) integer :: i,j,tmp_k,tmp_i, tmp_j double precision :: max_elem, t1,t2,t3 print*,'' print*,'---hessian_FB---' call wall_time(t1) ! Allocation allocate(beta(tmp_list_size,tmp_list_size)) ! Calculation do tmp_j = 1, tmp_list_size j = tmp_list(tmp_j) do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) beta(tmp_i,tmp_j) = (mo_dipole_x(i,i) - mo_dipole_x(j,j))**2 - 4d0 * mo_dipole_x(i,j)**2 & +(mo_dipole_y(i,i) - mo_dipole_y(j,j))**2 - 4d0 * mo_dipole_y(i,j)**2 & +(mo_dipole_z(i,i) - mo_dipole_z(j,j))**2 - 4d0 * mo_dipole_z(i,j)**2 enddo enddo ! Diagonal of the hessian H = 0d0 do tmp_k = 1, tmp_n call vec_to_mat_index(tmp_k,tmp_i,tmp_j) H(tmp_k) = 4d0 * beta(tmp_i, tmp_j) enddo ! Deallocation deallocate(beta) call wall_time(t2) t3 = t2 - t1 print*,'Time in hessian_FB:', t3 print*,'---End hessian_FB---' end subroutine #+END_SRC *** Hessian (OMP) #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine hessian_FB_omp(tmp_n, tmp_list_size, tmp_list, H) implicit none BEGIN_DOC ! Compute the diagonal hessian for the Foster-Boys localization END_DOC integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) double precision, intent(out) :: H(tmp_n) double precision, allocatable :: beta(:,:) integer :: i,j,tmp_k,tmp_i,tmp_j double precision :: max_elem, t1,t2,t3 print*,'' print*,'---hessian_FB_omp---' call wall_time(t1) ! Allocation allocate(beta(tmp_list_size,tmp_list_size)) ! Initialization omp call omp_set_max_active_levels(1) !$OMP PARALLEL & !$OMP PRIVATE(i,j,tmp_i,tmp_j,tmp_k) & !$OMP SHARED(tmp_n,tmp_list_size,beta,H,mo_dipole_x,mo_dipole_y,mo_dipole_z,tmp_list) & !$OMP DEFAULT(NONE) ! Calculation !$OMP DO do tmp_j = 1, tmp_list_size j = tmp_list(tmp_j) do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) beta(tmp_i,tmp_j) = (mo_dipole_x(i,i) - mo_dipole_x(j,j))**2 - 4d0 * mo_dipole_x(i,j)**2 & +(mo_dipole_y(i,i) - mo_dipole_y(j,j))**2 - 4d0 * mo_dipole_y(i,j)**2 & +(mo_dipole_z(i,i) - mo_dipole_z(j,j))**2 - 4d0 * mo_dipole_z(i,j)**2 enddo enddo !$OMP END DO ! Initialization !$OMP DO do i = 1, tmp_n H(i) = 0d0 enddo !$OMP END DO ! Diagonalm of the hessian !$OMP DO do tmp_k = 1, tmp_n call vec_to_mat_index(tmp_k,tmp_i,tmp_j) H(tmp_k) = 4d0 * beta(tmp_i, tmp_j) enddo !$OMP END DO !$OMP END PARALLEL call omp_set_max_active_levels(4) ! Deallocation deallocate(beta) call wall_time(t2) t3 = t2 - t1 print*,'Time in hessian_FB_omp:', t3 print*,'---End hessian_FB_omp---' end subroutine #+END_SRC ** Pipek-Mezey *** Gradient v1 #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine grad_pipek(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) implicit none BEGIN_DOC ! Compute gradient for the Pipek-Mezey localization END_DOC 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, allocatable :: m_grad(:,:), tmp_int(:,:) integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho ! Allocation allocate(m_grad(tmp_list_size, tmp_list_size), tmp_int(tmp_list_size, tmp_list_size)) ! Initialization m_grad = 0d0 do a = 1, nucl_num ! loop over the nuclei tmp_int = 0d0 ! Initialization for each nuclei ! Loop over the MOs of the a given mo_class to compute do tmp_j = 1, tmp_list_size j = tmp_list(tmp_j) do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) 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 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) & + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j)) enddo enddo enddo enddo ! Gradient do tmp_j = 1, tmp_list_size do tmp_i = 1, tmp_list_size m_grad(tmp_i,tmp_j) = m_grad(tmp_i,tmp_j) + 4d0 * tmp_int(tmp_i,tmp_j) * (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j)) enddo enddo enddo ! 2D -> 1D do tmp_k = 1, tmp_n call vec_to_mat_index(tmp_k,tmp_i,tmp_j) v_grad(tmp_k) = m_grad(tmp_i,tmp_j) enddo ! Maximum element in the gradient max_elem = 0d0 do tmp_k = 1, tmp_n if (ABS(v_grad(tmp_k)) > max_elem) then max_elem = ABS(v_grad(tmp_k)) endif enddo ! Norm of the gradient norm_grad = 0d0 do tmp_k = 1, tmp_n norm_grad = norm_grad + v_grad(tmp_k)**2 enddo norm_grad = dsqrt(norm_grad) print*, 'Maximal element in the gradient:', max_elem print*, 'Norm of the gradient:', norm_grad ! Deallocation deallocate(m_grad,tmp_int) end subroutine grad_pipek #+END_SRC *** Gradient The gradient is \begin{align*} \left. \frac{\partial \mathcal{P} (\theta)}{\partial \theta} \right|_{\theta=0}= \gamma^{PM} \end{align*} with \begin{align*} \gamma_{st}^{PM} = \sum_{A=1}^N \left[ - \right] \end{align*} \begin{align*} = \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*} $\sum_{\rho}$ -> sum over all the AOs $\sum_{\mu \in A}$ -> sum over the AOs which belongs to atom A $c^t$ -> expansion coefficient of orbital |t> Input: | 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(tmp_list_size) | integer | MOs in the mo_class | Output: | v_grad(tmp_n) | double precision | Gradient in the subspace | | max_elem | double precision | Maximal element in the gradient | | norm_grad | double precision | Norm of the gradient | Internal: | m_grad(tmp_list_size,tmp_list_size) | double precision | Gradient in a 2D array | | tmp_int(tmp_list_size,tmp_list_size) | | Temporary array to store the integrals | | tmp_accu(tmp_list_size,tmp_list_size) | | Temporary array to store a matrix | | | | product and compute tmp_int | | CS(tmp_list_size,ao_num) | | Array to store the result of mo_coef * ao_overlap | | tmp_mo_coef(ao_num,tmp_list_size) | | Array to store just the useful MO coefficients | | | | depending of the mo_class | | tmp_mo_coef2(nucl_n_aos(a),tmp_list_size) | | Array to store just the useful MO coefficients | | | | depending of the nuclei | | tmp_CS(tmp_list_size,nucl_n_aos(a)) | | Array to store just the useful mo_coef * ao_overlap | | | | values depending of the nuclei | | a | | index to loop over the nuclei | | b | | index to loop over the AOs which belongs to the nuclei a | | mu | | index to refer to an AO which belongs to the nuclei a | | rho | | index to loop over all the AOs | #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine gradient_PM(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad) implicit none BEGIN_DOC ! Compute gradient for the Pipek-Mezey localization END_DOC 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, allocatable :: m_grad(:,:), tmp_int(:,:), CS(:,:), tmp_mo_coef(:,:), tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:) integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho double precision :: t1,t2,t3 print*,'' print*,'---gradient_PM---' call wall_time(t1) ! Allocation allocate(m_grad(tmp_list_size, tmp_list_size), tmp_int(tmp_list_size, tmp_list_size),tmp_accu(tmp_list_size, tmp_list_size)) allocate(CS(tmp_list_size,ao_num),tmp_mo_coef(ao_num,tmp_list_size)) ! submatrix of the mo_coef do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) do j = 1, ao_num tmp_mo_coef(j,tmp_i) = mo_coef(j,i) 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)) m_grad = 0d0 do a = 1, nucl_num ! loop over the nuclei tmp_int = 0d0 !do tmp_j = 1, tmp_list_size ! do tmp_i = 1, tmp_list_size ! do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a ! mu = nucl_aos(a,b) ! tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (CS(tmp_i,mu) * tmp_mo_coef(mu,tmp_j) + tmp_mo_coef(mu,tmp_i) * CS(tmp_j,mu)) ! ! (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) & ! !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j)) ! enddo ! enddo !enddo 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 b = 1, nucl_n_aos(a) mu = nucl_aos(a,b) tmp_mo_coef2(b,tmp_i) = tmp_mo_coef(mu,tmp_i) enddo enddo do b = 1, nucl_n_aos(a) mu = nucl_aos(a,b) do tmp_i = 1, tmp_list_size tmp_CS(tmp_i,b) = CS(tmp_i,mu) 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)) do tmp_j = 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)) enddo enddo deallocate(tmp_mo_coef2,tmp_CS) do tmp_j = 1, tmp_list_size do tmp_i = 1, tmp_list_size m_grad(tmp_i,tmp_j) = m_grad(tmp_i,tmp_j) + 4d0 * tmp_int(tmp_i,tmp_j) * (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j)) enddo enddo enddo ! 2D -> 1D do tmp_k = 1, tmp_n call vec_to_mat_index(tmp_k,tmp_i,tmp_j) v_grad(tmp_k) = m_grad(tmp_i,tmp_j) enddo ! Maximum element in the gradient max_elem = 0d0 do tmp_k = 1, tmp_n if (ABS(v_grad(tmp_k)) > max_elem) then max_elem = ABS(v_grad(tmp_k)) endif enddo ! Norm of the gradient norm_grad = 0d0 do tmp_k = 1, tmp_n norm_grad = norm_grad + v_grad(tmp_k)**2 enddo norm_grad = dsqrt(norm_grad) print*, 'Maximal element in the gradient:', max_elem print*, 'Norm of the gradient:', norm_grad ! Deallocation deallocate(m_grad,tmp_int,CS,tmp_mo_coef) call wall_time(t2) t3 = t2 - t1 print*,'Time in gradient_PM:', t3 print*,'---End gradient_PM---' end #+END_SRC *** Hessian v1 #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine hess_pipek(tmp_n, tmp_list_size, tmp_list, H) implicit none BEGIN_DOC ! Compute diagonal hessian for the Pipek-Mezey localization END_DOC integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) double precision, intent(out) :: H(tmp_n) double precision, allocatable :: beta(:,:),tmp_int(:,:) integer :: i,j,tmp_k,tmp_i, tmp_j, a,b,rho,mu double precision :: max_elem ! Allocation allocate(beta(tmp_list_size,tmp_list_size),tmp_int(tmp_list_size,tmp_list_size)) beta = 0d0 do a = 1, nucl_num tmp_int = 0d0 do tmp_j = 1, tmp_list_size j = tmp_list(tmp_j) do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) do rho = 1, ao_num do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a mu = nucl_aos(a,b) 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)) enddo enddo enddo enddo ! Calculation do tmp_j = 1, tmp_list_size do tmp_i = 1, tmp_list_size beta(tmp_i,tmp_j) = beta(tmp_i, tmp_j) + (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j))**2 - 4d0 * tmp_int(tmp_i,tmp_j)**2 enddo enddo enddo H = 0d0 do tmp_k = 1, tmp_n call vec_to_mat_index(tmp_k,tmp_i,tmp_j) H(tmp_k) = 4d0 * beta(tmp_i, tmp_j) enddo ! Deallocation deallocate(beta,tmp_int) end #+END_SRC *** Hessian The hessian is \begin{align*} \left. \frac{\partial^2 \mathcal{P} (\theta)}{\partial \theta^2}\right|_{\theta=0} = 4 \beta^{PM} \end{align*} \begin{align*} \beta_{st}^{PM} = \sum_{A=1}^N \left( ^2 - \frac{1}{4} \left[ - \right]^2 \right) \end{align*} with \begin{align*} = \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*} $\sum_{\rho}$ -> sum over all the AOs $\sum_{\mu \in A}$ -> sum over the AOs which belongs to atom A $c^t$ -> expansion coefficient of orbital |t> #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H) implicit none BEGIN_DOC ! Compute diagonal hessian for the Pipek-Mezey localization END_DOC integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) double precision, intent(out) :: H(tmp_n) 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 double precision :: max_elem, t1,t2,t3 print*,'' print*,'---hessian_PM---' call wall_time(t1) ! Allocation allocate(beta(tmp_list_size,tmp_list_size),tmp_int(tmp_list_size,tmp_list_size),tmp_accu(tmp_list_size,tmp_list_size)) allocate(CS(tmp_list_size,ao_num),tmp_mo_coef(ao_num,tmp_list_size)) beta = 0d0 do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) do j = 1, ao_num tmp_mo_coef(j,tmp_i) = mo_coef(j,i) 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)) do a = 1, nucl_num ! loop over the nuclei tmp_int = 0d0 !do tmp_j = 1, tmp_list_size ! do tmp_i = 1, tmp_list_size ! do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a ! mu = nucl_aos(a,b) ! tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (CS(tmp_i,mu) * tmp_mo_coef(mu,tmp_j) + tmp_mo_coef(mu,tmp_i) * CS(tmp_j,mu)) ! ! (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) & ! !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j)) ! enddo ! enddo !enddo 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 b = 1, nucl_n_aos(a) mu = nucl_aos(a,b) tmp_mo_coef2(b,tmp_i) = tmp_mo_coef(mu,tmp_i) enddo enddo do b = 1, nucl_n_aos(a) mu = nucl_aos(a,b) do tmp_i = 1, tmp_list_size tmp_CS(tmp_i,b) = CS(tmp_i,mu) 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)) do tmp_j = 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)) enddo enddo deallocate(tmp_mo_coef2,tmp_CS) ! Calculation do tmp_j = 1, tmp_list_size do tmp_i = 1, tmp_list_size beta(tmp_i,tmp_j) = beta(tmp_i, tmp_j) + (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j))**2 - 4d0 * tmp_int(tmp_i,tmp_j)**2 enddo enddo enddo H = 0d0 do tmp_k = 1, tmp_n call vec_to_mat_index(tmp_k,tmp_i,tmp_j) H(tmp_k) = 4d0 * beta(tmp_i, tmp_j) enddo ! Deallocation deallocate(beta,tmp_int) call wall_time(t2) t3 = t2 - t1 print*,'Time in hessian_PM:', t3 print*,'---End hessian_PM---' end #+END_SRC ** Criterion *** Criterion PM (old) #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine compute_crit_pipek(criterion) implicit none BEGIN_DOC ! Compute the Pipek-Mezey localization criterion END_DOC double precision, intent(out) :: criterion double precision, allocatable :: tmp_int(:,:) integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho ! Allocation allocate(tmp_int(mo_num, mo_num)) criterion = 0d0 do a = 1, nucl_num ! loop over the nuclei tmp_int = 0d0 do i = 1, mo_num 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 mu = nucl_aos(a,b) tmp_int(i,i) = tmp_int(i,i) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,i) & + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,i)) enddo enddo enddo do i = 1, mo_num criterion = criterion + tmp_int(i,i)**2 enddo enddo criterion = - criterion deallocate(tmp_int) end #+END_SRC *** Criterion PM The criterion is computed as \begin{align*} \mathcal{P} = \sum_{i=1}^n \sum_{A=1}^N \left[ \right]^2 \end{align*} with \begin{align*} = \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*} #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine criterion_PM(tmp_list_size,tmp_list,criterion) implicit none BEGIN_DOC ! Compute the Pipek-Mezey localization criterion END_DOC integer, intent(in) :: tmp_list_size, tmp_list(tmp_list_size) double precision, intent(out) :: criterion double precision, allocatable :: tmp_int(:,:),CS(:,:) integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho print*,'' print*,'---criterion_PM---' ! Allocation allocate(tmp_int(tmp_list_size, tmp_list_size),CS(mo_num,ao_num)) ! Initialization 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)) do a = 1, nucl_num ! loop over the nuclei tmp_int = 0d0 do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a 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)) ! (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) & !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j)) enddo enddo do tmp_i = 1, tmp_list_size criterion = criterion + tmp_int(tmp_i,tmp_i)**2 enddo enddo criterion = - criterion deallocate(tmp_int,CS) print*,'---End criterion_PM---' end #+END_SRC *** Criterion PM v3 #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion) implicit none BEGIN_DOC ! Compute the Pipek-Mezey localization criterion END_DOC integer, intent(in) :: tmp_list_size, tmp_list(tmp_list_size) double precision, intent(out) :: criterion double precision, allocatable :: tmp_int(:,:), CS(:,:), tmp_mo_coef(:,:), tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:) integer :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho,nu,c double precision :: t1,t2,t3 print*,'' print*,'---criterion_PM_v3---' call wall_time(t1) ! Allocation allocate(tmp_int(tmp_list_size, tmp_list_size),tmp_accu(tmp_list_size, tmp_list_size)) allocate(CS(tmp_list_size,ao_num),tmp_mo_coef(ao_num,tmp_list_size)) criterion = 0d0 ! submatrix of the mo_coef do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) do j = 1, ao_num tmp_mo_coef(j,tmp_i) = mo_coef(j,i) enddo enddo ! ao_overlap(ao_num,ao_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)) do a = 1, nucl_num ! loop over the nuclei do j = 1, tmp_list_size do i = 1, tmp_list_size tmp_int(i,j) = 0d0 enddo enddo !do tmp_j = 1, tmp_list_size ! do tmp_i = 1, tmp_list_size ! do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a ! mu = nucl_aos(a,b) ! tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (CS(tmp_i,mu) * tmp_mo_coef(mu,tmp_j) + tmp_mo_coef(mu,tmp_i) * CS(tmp_j,mu)) ! ! (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) & ! !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j)) ! enddo ! enddo !enddo 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 b = 1, nucl_n_aos(a) mu = nucl_aos(a,b) tmp_mo_coef2(b,tmp_i) = tmp_mo_coef(mu,tmp_i) enddo enddo do b = 1, nucl_n_aos(a) mu = nucl_aos(a,b) do tmp_i = 1, tmp_list_size tmp_CS(tmp_i,b) = CS(tmp_i,mu) 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)) ! Integrals do tmp_j = 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)) enddo enddo deallocate(tmp_mo_coef2,tmp_CS) ! Criterion do tmp_i = 1, tmp_list_size criterion = criterion + tmp_int(tmp_i,tmp_i)**2 enddo enddo criterion = - criterion deallocate(tmp_int,CS,tmp_accu,tmp_mo_coef) call wall_time(t2) t3 = t2 - t1 print*,'Time in criterion_PM_v3:', t3 print*,'---End criterion_PM_v3---' end #+END_SRC *** Criterion FB (old) The criterion is just computed as \begin{align*} C = - \sum_i^{mo_{num}} (^2 + ^2 + ^2) \end{align*} The minus sign is here in order to minimize this criterion Output: | criterion | double precision | criterion for the Foster-Boys localization | #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine criterion_FB_old(criterion) implicit none BEGIN_DOC ! Compute the Foster-Boys localization criterion END_DOC double precision, intent(out) :: criterion integer :: i ! Criterion (= \sum_i ^2 ) criterion = 0d0 do i = 1, mo_num criterion = criterion + mo_dipole_x(i,i)**2 + mo_dipole_y(i,i)**2 + mo_dipole_z(i,i)**2 enddo criterion = - criterion end subroutine #+END_SRC *** Criterion FB #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine criterion_FB(tmp_list_size, tmp_list, criterion) implicit none BEGIN_DOC ! Compute the Foster-Boys localization criterion END_DOC integer, intent(in) :: tmp_list_size, tmp_list(tmp_list_size) double precision, intent(out) :: criterion integer :: i, tmp_i ! Criterion (= - \sum_i ^2 ) criterion = 0d0 do tmp_i = 1, tmp_list_size i = tmp_list(tmp_i) criterion = criterion + mo_dipole_x(i,i)**2 + mo_dipole_y(i,i)**2 + mo_dipole_z(i,i)**2 enddo criterion = - criterion end subroutine #+END_SRC ** Theta In: | n | integer | number of MOs in the considered MO class | | l | integer | list of MOs of the considered class | Out: | m_x(n,n) | double precision | Matrix containing the rotation angle between all the different | | | | pairs of MOs to apply the rotations (need a minus sign) | | max_elem | double precision | Maximal angle in absolute value | $$\cos(4 \theta) = \frac{-A{ij}}{\sqrt{(A_{ij}^2 + B_{ij}^2)} $$ $$\sin(4 \theta) = \frac{B{ij}}{\sqrt{(A_{ij}^2 + B_{ij}^2)} $$ $$\tan(4 \theta) = \frac{\sin(4 \theta)}{\cos(4 \theta)}$$ where $\theta$ is in fact $\theta_{ij}$ For Foster-Boys localization: $$A_{ij} = ^2 - \frac{1}{4} ( - )^2$$ $$B_{ij} = ( - )$$ For Pipek-Mezey localization: $$A_{ij} = \sum_A ^2 - \frac{1}{4} ( - )^2$$ $$B_{ij} = \sum_A ( - )$$ with $$ = \frac{1}{2} \sum_\rho \sum_{\mu \in A} ( c_\rho^{i*} S_{\rho \mu} c_\mu^j + c_\mu^{i*} S_{\mu \rho} c_\rho^j)$$ $i,j$ MOs $\mu, \rho$ AOs $A$ nucleus $S$ overlap matrix $c$ MO coefficient $r$ position operator #+begin_src f90 :tangle localization_sub.irp.f subroutine theta_FB(l, n, m_x, max_elem) include 'pi.h' BEGIN_DOC ! Compute the angles to minimize the Foster-Boys criterion by using pairwise rotations of the MOs ! Warning: you must give - the angles to build the rotation matrix... END_DOC implicit none integer, intent(in) :: n, l(n) double precision, intent(out) :: m_x(n,n), max_elem integer :: i,j, tmp_i, tmp_j double precision, allocatable :: cos4theta(:,:), sin4theta(:,:) double precision, allocatable :: A(:,:), B(:,:), beta(:,:), gamma(:,:) integer :: idx_i,idx_j allocate(cos4theta(n, n), sin4theta(n, n)) allocate(A(n,n), B(n,n), beta(n,n), gamma(n,n)) do tmp_j = 1, n j = l(tmp_j) do tmp_i = 1, n i = l(tmp_i) A(tmp_i,tmp_j) = mo_dipole_x(i,j)**2 - 0.25d0 * (mo_dipole_x(i,i) - mo_dipole_x(j,j))**2 & + mo_dipole_y(i,j)**2 - 0.25d0 * (mo_dipole_y(i,i) - mo_dipole_y(j,j))**2 & + mo_dipole_z(i,j)**2 - 0.25d0 * (mo_dipole_z(i,i) - mo_dipole_z(j,j))**2 enddo A(j,j) = 0d0 enddo do tmp_j = 1, n j = l(tmp_j) do tmp_i = 1, n i = l(tmp_i) B(tmp_i,tmp_j) = mo_dipole_x(i,j) * (mo_dipole_x(i,i) - mo_dipole_x(j,j)) & + mo_dipole_y(i,j) * (mo_dipole_y(i,i) - mo_dipole_y(j,j)) & + mo_dipole_z(i,j) * (mo_dipole_z(i,i) - mo_dipole_z(j,j)) enddo enddo !do tmp_j = 1, n ! j = l(tmp_j) ! do tmp_i = 1, n ! i = l(tmp_i) ! beta(tmp_i,tmp_j) = (mo_dipole_x(i,i) - mo_dipole_x(j,j)) - 4d0 * mo_dipole_x(i,j)**2 & ! + (mo_dipole_y(i,i) - mo_dipole_y(j,j)) - 4d0 * mo_dipole_y(i,j)**2 & ! + (mo_dipole_z(i,i) - mo_dipole_z(j,j)) - 4d0 * mo_dipole_z(i,j)**2 ! enddo !enddo !do tmp_j = 1, n ! j = l(tmp_j) ! do tmp_i = 1, n ! i = l(tmp_i) ! gamma(tmp_i,tmp_j) = 4d0 * ( mo_dipole_x(i,j) * (mo_dipole_x(i,i) - mo_dipole_x(j,j)) & ! + mo_dipole_y(i,j) * (mo_dipole_y(i,i) - mo_dipole_y(j,j)) & ! + mo_dipole_z(i,j) * (mo_dipole_z(i,i) - mo_dipole_z(j,j))) ! enddo !enddo ! !do j = 1, n ! do i = 1, n ! cos4theta(i,j) = - A(i,j) / dsqrt(A(i,j)**2 + B(i,j)**2) ! enddo !enddo !do j = 1, n ! do i = 1, n ! sin4theta(i,j) = B(i,j) / dsqrt(A(i,j)**2 + B(i,j)**2) ! enddo !enddo ! Theta do j = 1, n do i = 1, n m_x(i,j) = 0.25d0 * atan2(B(i,j), -A(i,j)) !m_x(i,j) = 0.25d0 * atan2(sin4theta(i,j), cos4theta(i,j)) enddo enddo ! Enforce a perfect antisymmetry do j = 1, n-1 do i = j+1, n m_x(j,i) = - m_x(i,j) enddo enddo do i = 1, n m_x(i,i) = 0d0 enddo ! Max max_elem = 0d0 do j = 1, n-1 do i = j+1, n if (dabs(m_x(i,j)) > dabs(max_elem)) then max_elem = m_x(i,j) !idx_i = i !idx_j = j endif enddo enddo ! Debug !print*,'' !print*,'sin/B' !do i = 1, n ! write(*,'(100F10.4)') sin4theta(i,:) ! !B(i,:) !enddo !print*,'cos/A' !do i = 1, n ! write(*,'(100F10.4)') cos4theta(i,:) ! !A(i,:) !enddo !print*,'X' !!m_x = 0d0 !!m_x(idx_i,idx_j) = max_elem !!m_x(idx_j,idx_i) = -max_elem !do i = 1, n ! write(*,'(100F10.4)') m_x(i,:) !enddo !print*,idx_i,idx_j,max_elem max_elem = dabs(max_elem) deallocate(cos4theta, sin4theta) deallocate(A,B,beta,gamma) end #+end_src #+begin_src f90 :comments org :tangle localization_sub.irp.f subroutine theta_PM(l, n, m_x, max_elem) include 'pi.h' BEGIN_DOC ! Compute the angles to minimize the Foster-Boys criterion by using pairwise rotations of the MOs ! Warning: you must give - the angles to build the rotation matrix... END_DOC implicit none integer, intent(in) :: n, l(n) double precision, intent(out) :: m_x(n,n), max_elem integer :: a,b,i,j,tmp_i,tmp_j,rho,mu,nu,idx_i,idx_j double precision, allocatable :: Aij(:,:), Bij(:,:), Pa(:,:) allocate(Aij(n,n), Bij(n,n), Pa(n,n)) do a = 1, nucl_num ! loop over the nuclei Pa = 0d0 ! Initialization for each nuclei ! Loop over the MOs of the a given mo_class to compute do tmp_j = 1, n j = l(tmp_j) do tmp_i = 1, n i = l(tmp_i) 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 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) & + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j)) enddo enddo enddo enddo ! A do j = 1, n do i = 1, n Aij(i,j) = Aij(i,j) + Pa(i,j)**2 - 0.25d0 * (Pa(i,i) - Pa(j,j))**2 enddo enddo ! B do j = 1, n do i = 1, n Bij(i,j) = Bij(i,j) + Pa(i,j) * (Pa(i,i) - Pa(j,j)) enddo enddo enddo ! Theta do j = 1, n do i = 1, n m_x(i,j) = 0.25d0 * atan2(Bij(i,j), -Aij(i,j)) enddo enddo ! Enforce a perfect antisymmetry do j = 1, n-1 do i = j+1, n m_x(j,i) = - m_x(i,j) enddo enddo do i = 1, n m_x(i,i) = 0d0 enddo ! Max max_elem = 0d0 do j = 1, n-1 do i = j+1, n if (dabs(m_x(i,j)) > dabs(max_elem)) then max_elem = m_x(i,j) idx_i = i idx_j = j endif enddo enddo ! Debug !do i = 1, n ! write(*,'(100F10.4)') m_x(i,:) !enddo !print*,'Max',idx_i,idx_j,max_elem max_elem = dabs(max_elem) deallocate(Aij,Bij,Pa) end #+end_src ** Spatial extent The spatial extent of an orbital $i$ is computed as \begin{align*} \sum_{\lambda=x,y,z}\sqrt{ - ^2} \end{align*} From that we can also compute the average and the standard deviation #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine compute_spatial_extent(spatial_extent) implicit none BEGIN_DOC ! Compute the spatial extent of the MOs END_DOC double precision, intent(out) :: spatial_extent(mo_num) double precision :: average_core, average_act, average_inact, average_virt double precision :: std_var_core, std_var_act, std_var_inact, std_var_virt integer :: i,j,k,l spatial_extent = 0d0 do i = 1, mo_num spatial_extent(i) = mo_spread_x(i,i) - mo_dipole_x(i,i)**2 enddo do i = 1, mo_num spatial_extent(i) = spatial_extent(i) + mo_spread_y(i,i) - mo_dipole_y(i,i)**2 enddo do i = 1, mo_num spatial_extent(i) = spatial_extent(i) + mo_spread_z(i,i) - mo_dipole_z(i,i)**2 enddo do i = 1, mo_num spatial_extent(i) = dsqrt(spatial_extent(i)) enddo average_core = 0d0 std_var_core = 0d0 if (dim_list_core_orb >= 2) then call compute_average_sp_ext(spatial_extent, list_core, dim_list_core_orb, average_core) call compute_std_var_sp_ext(spatial_extent, list_core, dim_list_core_orb, average_core, std_var_core) endif average_act = 0d0 std_var_act = 0d0 if (dim_list_act_orb >= 2) then 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) endif average_inact = 0d0 std_var_inact = 0d0 if (dim_list_inact_orb >= 2) then call compute_average_sp_ext(spatial_extent, list_inact, dim_list_inact_orb, average_inact) call compute_std_var_sp_ext(spatial_extent, list_inact, dim_list_inact_orb, average_inact, std_var_inact) endif average_virt = 0d0 std_var_virt = 0d0 if (dim_list_virt_orb >= 2) then call compute_average_sp_ext(spatial_extent, list_virt, dim_list_virt_orb, average_virt) call compute_std_var_sp_ext(spatial_extent, list_virt, dim_list_virt_orb, average_virt, std_var_virt) endif print*,'' print*,'=============================' print*,' Spatial extent of the MOs' print*,'=============================' print*,'' print*, 'elec_num:', elec_num print*, 'elec_alpha_num:', elec_alpha_num print*, 'elec_beta_num:', elec_beta_num print*, 'core:', dim_list_core_orb print*, 'act:', dim_list_act_orb print*, 'inact:', dim_list_inact_orb print*, 'virt:', dim_list_virt_orb print*, 'mo_num:', mo_num print*,'' print*,'-- Core MOs --' print*,'Average:', average_core print*,'Std var:', std_var_core print*,'' print*,'-- Active MOs --' print*,'Average:', average_act print*,'Std var:', std_var_act print*,'' print*,'-- Inactive MOs --' print*,'Average:', average_inact print*,'Std var:', std_var_inact print*,'' print*,'-- Virtual MOs --' print*,'Average:', average_virt print*,'Std var:', std_var_virt print*,'' print*,'Spatial extent:' do i = 1, mo_num print*, i, spatial_extent(i) enddo end #+END_SRC #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine compute_average_sp_ext(spatial_extent, list, list_size, average) implicit none BEGIN_DOC ! Compute the average spatial extent of the MOs END_DOC integer, intent(in) :: list_size, list(list_size) double precision, intent(in) :: spatial_extent(mo_num) double precision, intent(out) :: average integer :: i, tmp_i average = 0d0 do tmp_i = 1, list_size i = list(tmp_i) average = average + spatial_extent(i) enddo average = average / DBLE(list_size) end #+END_SRC #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine compute_std_var_sp_ext(spatial_extent, list, list_size, average, std_var) implicit none BEGIN_DOC ! Compute the standard deviation of the spatial extent of the MOs END_DOC integer, intent(in) :: list_size, list(list_size) double precision, intent(in) :: spatial_extent(mo_num) double precision, intent(in) :: average double precision, intent(out) :: std_var integer :: i, tmp_i std_var = 0d0 do tmp_i = 1, list_size i = list(tmp_i) std_var = std_var + (spatial_extent(i) - average)**2 enddo std_var = dsqrt(1d0/DBLE(list_size) * std_var) end #+END_SRC ** Utils #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine apply_pre_rotation() implicit none BEGIN_DOC ! Apply a rotation between the MOs END_DOC double precision, allocatable :: pre_rot(:,:), prev_mos(:,:), R(:,:) double precision :: t1,t2,t3 integer :: i,j,tmp_i,tmp_j integer :: info logical :: enforce_step_cancellation print*,'---apply_pre_rotation---' call wall_time(t1) allocate(pre_rot(mo_num,mo_num), prev_mos(ao_num,mo_num), R(mo_num,mo_num)) ! Initialization of the matrix pre_rot = 0d0 if (kick_in_mos) then ! Pre rotation for core MOs if (dim_list_core_orb >= 2) then do tmp_j = 1, dim_list_core_orb j = list_core(tmp_j) do tmp_i = 1, dim_list_core_orb i = list_core(tmp_i) if (i > j) then pre_rot(i,j) = angle_pre_rot elseif (i < j) then pre_rot(i,j) = - angle_pre_rot else pre_rot(i,j) = 0d0 endif enddo enddo endif ! Pre rotation for active MOs if (dim_list_act_orb >= 2) then do tmp_j = 1, dim_list_act_orb j = list_act(tmp_j) do tmp_i = 1, dim_list_act_orb i = list_act(tmp_i) if (i > j) then pre_rot(i,j) = angle_pre_rot elseif (i < j) then pre_rot(i,j) = - angle_pre_rot else pre_rot(i,j) = 0d0 endif enddo enddo endif ! Pre rotation for inactive MOs if (dim_list_inact_orb >= 2) then do tmp_j = 1, dim_list_inact_orb j = list_inact(tmp_j) do tmp_i = 1, dim_list_inact_orb i = list_inact(tmp_i) if (i > j) then pre_rot(i,j) = angle_pre_rot elseif (i < j) then pre_rot(i,j) = - angle_pre_rot else pre_rot(i,j) = 0d0 endif enddo enddo endif ! Pre rotation for virtual MOs if (dim_list_virt_orb >= 2) then do tmp_j = 1, dim_list_virt_orb j = list_virt(tmp_j) do tmp_i = 1, dim_list_virt_orb i = list_virt(tmp_i) if (i > j) then pre_rot(i,j) = angle_pre_rot elseif (i < j) then pre_rot(i,j) = - angle_pre_rot else pre_rot(i,j) = 0d0 endif enddo enddo endif ! Nothing for deleted ones ! Compute pre rotation matrix from pre_rot call rotation_matrix(pre_rot,mo_num,R,mo_num,mo_num,info,enforce_step_cancellation) if (enforce_step_cancellation) then print*, 'Cancellation of the pre rotation, too big error in the rotation matrix' print*, 'Reduce the angle for the pre rotation, abort' call abort endif ! New Mos (we don't car eabout the previous MOs prev_mos) call apply_mo_rotation(R,prev_mos) ! Update the things related to mo_coef TOUCH mo_coef call save_mos endif deallocate(pre_rot, prev_mos, R) call wall_time(t2) t3 = t2-t1 print*,'Time in apply_pre_rotation:', t3 print*,'---End apply_pre_rotation---' end #+END_SRC #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine x_tmp_orb_loc_v2(tmp_n, tmp_list_size, tmp_list, v_grad, H,tmp_x, tmp_m_x) implicit none integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) double precision, intent(in) :: v_grad(tmp_n) double precision, intent(in) :: H(tmp_n, tmp_n) double precision, intent(out) :: tmp_m_x(tmp_list_size, tmp_list_size), tmp_x(tmp_list_size) !double precision, allocatable :: x(:) double precision :: lambda , accu, max_elem integer :: i,j,tmp_i,tmp_j,tmp_k ! Allocation !allocate(x(tmp_n)) ! Level shifted hessian lambda = 0d0 do tmp_k = 1, tmp_n if (H(tmp_k,tmp_k) < lambda) then lambda = H(tmp_k,tmp_k) endif enddo ! min element in the hessian if (lambda < 0d0) then lambda = -lambda + 1d-6 endif print*, 'lambda', lambda ! Good do tmp_k = 1, tmp_n 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)) !x(tmp_k) = - 1d0/(ABS(H(tmp_k,tmp_k))+lambda) * (-v_grad(tmp_k)) endif enddo ! 1D tmp -> 2D tmp tmp_m_x = 0d0 do tmp_j = 1, tmp_list_size - 1 do tmp_i = tmp_j + 1, tmp_list_size call mat_to_vec_index(tmp_i,tmp_j,tmp_k) tmp_m_x(tmp_i, tmp_j) = tmp_x(tmp_k)!x(tmp_k) enddo enddo ! Antisym do tmp_i = 1, tmp_list_size - 1 do tmp_j = tmp_i + 1, tmp_list_size tmp_m_x(tmp_i,tmp_j) = - tmp_m_x(tmp_j,tmp_i) enddo enddo ! Deallocation !deallocate(x) end subroutine #+END_SRC #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine ao_to_mo_no_sym(A_ao,LDA_ao,A_mo,LDA_mo) implicit none BEGIN_DOC ! Transform A from the |AO| basis to the |MO| basis ! ! $C^\dagger.A_{ao}.C$ END_DOC integer, intent(in) :: LDA_ao,LDA_mo double precision, intent(in) :: A_ao(LDA_ao,ao_num) double precision, intent(out) :: A_mo(LDA_mo,mo_num) double precision, allocatable :: T(:,:) allocate ( T(ao_num,mo_num) ) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T call dgemm('N','N', ao_num, mo_num, ao_num, & 1.d0, A_ao,LDA_ao, & mo_coef, size(mo_coef,1), & 0.d0, T, size(T,1)) call dgemm('T','N', mo_num, mo_num, ao_num, & 1.d0, mo_coef,size(mo_coef,1), & T, ao_num, & 0.d0, A_mo, size(A_mo,1)) deallocate(T) end #+END_SRC #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f subroutine run_sort_by_fock_energies() implicit none BEGIN_DOC ! Saves the current MOs ordered by diagonal element of the Fock operator. END_DOC integer :: i,j,k,l,tmp_i,tmp_k,tmp_list_size integer, allocatable :: iorder(:), tmp_list(:) double precision, allocatable :: fock_energies_tmp(:), tmp_mo_coef(:,:) ! Test do l = 1, 4 if (l==1) then ! core tmp_list_size = dim_list_core_orb elseif (l==2) then ! act tmp_list_size = dim_list_act_orb elseif (l==3) then ! inact tmp_list_size = dim_list_inact_orb else ! virt tmp_list_size = dim_list_virt_orb endif if (tmp_list_size >= 2) then ! Allocation tmp array allocate(tmp_list(tmp_list_size)) ! To give the list of MOs in a mo_class if (l==1) then ! core tmp_list = list_core elseif (l==2) then tmp_list = list_act elseif (l==3) then tmp_list = list_inact else tmp_list = list_virt endif 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)) !print*,'MOs before sorting them by f_p^p energies:' do i = 1, tmp_list_size tmp_i = tmp_list(i) fock_energies_tmp(i) = Fock_matrix_diag_mo(tmp_i) iorder(i) = i !print*, tmp_i, fock_energies_tmp(i) enddo call dsort(fock_energies_tmp, iorder, tmp_list_size) print*,'MOs after sorting them by f_p^p energies:' do i = 1, tmp_list_size k = iorder(i) tmp_k = tmp_list(k) print*, tmp_k, fock_energies_tmp(k) do j = 1, ao_num tmp_mo_coef(j,k) = mo_coef(j,tmp_k) enddo enddo ! Update the MOs after sorting them by energies do i = 1, tmp_list_size tmp_i = tmp_list(i) do j = 1, ao_num mo_coef(j,tmp_i) = tmp_mo_coef(j,i) enddo enddo if (debug_hf) then touch mo_coef print*,'HF energy:', HF_energy endif print*,'' deallocate(iorder, fock_energies_tmp, tmp_list, tmp_mo_coef) endif enddo touch mo_coef call save_mos end #+END_SRC #+BEGIN_SRC f90 :comments org :tangle localization_sub.irp.f function is_core(i) implicit none BEGIN_DOC ! True if the orbital i is a core orbital END_DOC integer, intent(in) :: i logical :: is_core integer :: j ! Init is_core = .False. ! Search do j = 1, dim_list_core_orb if (list_core(j) == i) then is_core = .True. exit endif enddo end function is_del(i) implicit none BEGIN_DOC ! True if the orbital i is a deleted orbital END_DOC integer, intent(in) :: i logical :: is_del integer :: j ! Init is_del = .False. ! Search do j = 1, dim_list_core_orb if (list_core(j) == i) then is_del = .True. exit endif enddo end subroutine set_classes_loc() implicit none integer :: i logical :: ok1, ok2 logical :: is_core, is_del integer(bit_kind) :: res(N_int,2) if (auto_mo_class) then do i = 1, mo_num if (is_core(i)) cycle if (is_del(i)) cycle call apply_hole(psi_det(1,1,1), 1, i, res, ok1, N_int) call apply_hole(psi_det(1,1,1), 2, i, res, ok2, N_int) if (ok1 .and. ok2) then mo_class(i) = 'Inactive' else if (.not. ok1 .and. .not. ok2) then mo_class(i) = 'Virtual' else mo_class(i) = 'Active' endif enddo touch mo_class endif end subroutine unset_classes_loc() implicit none integer :: i logical :: ok1, ok2 logical :: is_core, is_del integer(bit_kind) :: res(N_int,2) if (auto_mo_class) then do i = 1, mo_num if (is_core(i)) cycle if (is_del(i)) cycle mo_class(i) = 'Active' enddo touch mo_class endif end #+END_SRC