init commit mo_localization

src/mo_localization/NEED

@@ -0,0 +1,2 @@
+hartree_fock
+utils_trust_region

src/mo_localization/README.md

@@ -0,0 +1,94 @@
+Please, use the ifort compiler  
+  
+Some parameters can be changed with qp edit in the utils_trust_region section 
+ 
+If you modify the .org files, don't forget to do (you need emacs):  
+``` 
+./TANGLE_org_mode.sh  
+ninja  
+```  
+
+The documentation can be read using:  
+Ctrl-C Ctrl-e l p  
+after opening the filename.org in emacs. It will produce a  
+filename.pdf.  
+(Not available for all the files)  
+!!! Warning: the documentation can contain some errors !!! 
+
+# Orbital localisation
+To localize the MOs:  
+```
+qp run localization  
+```
+After that the ezfio directory contains the localized MOs  
+ 
+But the mo_class must be defined before, run 
+```
+qp set_mo_class -q
+```
+for more information or  
+```
+qp set_mo_class -c [] -a [] -v [] -i [] -d [] 
+```
+to set the mo classes. We don't care about the name of the   
+mo classes. The algorithm just localizes all the MOs of  
+a given class between them, for all the classes, except the deleted MOs.  
+
+If you just on kind of mo class to localize all the MOs between them  
+you have to put:
+```
+qp set mo_localization security_mo_class false
+```
+
+Before the localization, a kick is done for each mo class  
+(except the deleted ones) to break the MOs. This is done by   
+doing a given rotation between the MOs.
+This feature can be removed by setting:
+```
+qp set mo_localization kick_in_mos false
+```
+and the default angle for the rotation can be changed with:
+```
+qp set mo_localization angle_pre_rot 1e-3 # or something else
+```
+
+After the localization, the MOs of each class (except the deleted ones)  
+can be sorted between them using the diagonal elements of  
+the fock matrix with:
+```
+qp set mo_localization sort_mos_by_e true
+```
+
+You can check the Hartree-Fock energy before/during/after the localization  
+by putting (only for debugging):
+```
+qp set mo_localization debug_hf true 
+```
+
+## Foster-Boys & Pipek-Mezey
+Foster-Boys:  
+``` 
+qp set mo_localization localization_method boys 
+``` 
+ 
+Pipek-Mezey:  
+``` 
+qp set mo_localization localization_method pipek 
+``` 
+
+# Break the spatial symmetry of the MOs
+To break the spatial symmetry of the MOs:   
+```
+qp run break_spatial_sym
+```
+The default angle for the rotations is too big for this kind of
+application, a value between 1e-3 and 1e-6 should break the spatial
+symmetry with just a small change in the energy:
+```
+qp set mo_localization angle_pre_rot 1e-3
+``` 
+
+# Further improvements: 
+- Cleaner repo 
+- Correction of the errors in the documentations 
+- option with/without trust region 

src/mo_localization/TANGLE_org_mode.sh

@@ -0,0 +1,7 @@
+#!/bin/sh   
+
+list='ls *.org'
+for element in $list    
+do   
+		emacs --batch $element -f org-babel-tangle
+done

src/mo_localization/break_spatial_sym.org

@@ -0,0 +1,43 @@
+! A small program to break the spatial symmetry of the MOs.
+
+! You have to defined your MO classes or set security_mo_class to false
+! with:
+! qp set orbital_optimization security_mo_class false
+
+! The default angle for the rotations is too big for this kind of
+! application, a value between 1e-3 and 1e-6 should break the spatial
+! symmetry with just a small change in the energy. 
+
+#+BEGIN_SRC f90 :comments org :tangle break_spatial_sym.irp.f
+program break_spatial_sym
+
+  !BEGIN_DOC
+  ! Break the symmetry of the MOs with a rotation
+  !END_DOC
+
+  implicit none
+
+  kick_in_mos = .True.
+  TOUCH kick_in_mos
+
+  print*, 'Security mo_class:', security_mo_class
+
+  ! The default mo_classes are setted only if the MOs to localize are not specified
+  if (security_mo_class .and. (dim_list_act_orb == mo_num .or. &
+      dim_list_core_orb + dim_list_act_orb == mo_num)) then
+
+    print*, 'WARNING'
+    print*, 'You must set different mo_class with qp set_mo_class'
+    print*, 'If you want to kick all the orbitals:'
+    print*, 'qp set orbital_optimization security_mo_class false'
+    print*, ''
+    print*, 'abort'
+
+    call abort
+  
+  endif
+  
+  call apply_pre_rotation
+  
+end
+#+END_SRC

src/mo_localization/debug_gradient_loc.org

@@ -0,0 +1,64 @@
+#+BEGIN_SRC f90 :comments org :tangle debug_gradient_loc.irp.f
+program debug_gradient_loc
+
+  !BEGIN_DOC
+  ! Check if the gradient is correct
+  !END_DOC
+
+  implicit none
+
+  integer :: list_size, n
+  integer, allocatable :: list(:)
+  double precision, allocatable :: v_grad(:), v_grad2(:)
+  double precision :: norm, max_elem, threshold, max_error
+  integer :: i, nb_error
+
+  threshold = 1d-12
+
+  list = list_act
+  list_size = dim_list_act_orb
+
+  n = list_size*(list_size-1)/2
+  
+  allocate(v_grad(n),v_grad2(n))
+
+  if (localization_method == 'boys') then
+    print*,'Foster-Boys'
+    call gradient_FB(n,list_size,list,v_grad,max_elem,norm)
+    call gradient_FB_omp(n,list_size,list,v_grad2,max_elem,norm)
+  elseif (localization_method == 'pipek') then
+    print*,'Pipek-Mezey'
+    call gradient_PM(n,list_size,list,v_grad,max_elem,norm)
+    call gradient_PM(n,list_size,list,v_grad2,max_elem,norm) 
+  else
+    print*,'Unknown localization_method, please select boys or pipek'
+    call abort
+  endif
+ 
+  do i = 1, n
+    print*,i,v_grad(i)
+  enddo
+
+  v_grad = v_grad - v_grad2
+
+  nb_error = 0
+  max_elem = 0d0
+
+  do i = 1, n
+    if (dabs(v_grad(i)) > threshold) then
+      print*,v_grad(i)
+      nb_error = nb_error + 1
+      if (dabs(v_grad(i)) > max_elem) then
+        max_elem = v_grad(i)
+      endif
+    endif
+  enddo
+
+  print*,'Threshold error', threshold
+  print*, 'Nb error', nb_error
+  print*,'Max error', max_elem
+
+  deallocate(v_grad,v_grad2)
+ 
+end
+#+END_SRC

src/mo_localization/debug_hessian_loc.org

@@ -0,0 +1,64 @@
+#+BEGIN_SRC f90 :comments org :tangle debug_hessian_loc.irp.f
+program debug_hessian_loc
+
+  !BEGIN_DOC
+  ! Check if the hessian is correct
+  !END_DOC
+  
+  implicit none
+
+  integer :: list_size, n
+  integer, allocatable :: list(:)
+  double precision, allocatable :: H(:,:), H2(:,:)
+  double precision :: threshold, max_error, max_elem
+  integer :: i, nb_error
+
+  threshold = 1d-12
+
+  list = list_act
+  list_size = dim_list_act_orb
+
+  n = list_size*(list_size-1)/2
+  
+  allocate(H(n,n),H2(n,n))
+
+  if (localization_method == 'boys') then
+    print*,'Foster-Boys'
+    call hessian_FB(n,list_size,list,H)
+    call hessian_FB_omp(n,list_size,list,H2)
+  elseif(localization_method == 'pipek') then
+    print*,'Pipek-Mezey'
+    call hessian_PM(n,list_size,list,H)
+    call hessian_PM(n,list_size,list,H2)
+  else
+    print*,'Unknown localization_method, please select boys or pipek'
+    call abort
+  endif
+ 
+  do i = 1, n
+    print*,i,H(i,i)
+  enddo
+
+  H = H - H2
+
+  nb_error = 0
+  max_elem = 0d0
+
+  do i = 1, n
+    if (dabs(H(i,i)) > threshold) then
+      print*,H(i,i)
+      nb_error = nb_error + 1
+      if (dabs(H(i,i)) > max_elem) then
+        max_elem = H(i,i)
+      endif
+    endif
+  enddo
+
+  print*,'Threshold error', threshold
+  print*, 'Nb error', nb_error
+  print*,'Max error', max_elem
+
+  deallocate(H,H2)
+  
+end
+#+END_SRC

src/mo_localization/kick_the_mos.org

@@ -0,0 +1,33 @@
+#+BEGIN_SRC f90 :comments org :tangle kick_the_mos.irp.f
+program kick_the_mos
+
+  !BEGIN_DOC
+  ! To do a small rotation of the MOs
+  !END_DOC
+
+  implicit none
+
+  kick_in_mos = .True.
+  TOUCH kick_in_mos
+
+  print*, 'Security mo_class:', security_mo_class
+
+  ! The default mo_classes are setted only if the MOs to localize are not specified
+  if (security_mo_class .and. (dim_list_act_orb == mo_num .or. &
+      dim_list_core_orb + dim_list_act_orb == mo_num)) then
+
+    print*, 'WARNING'
+    print*, 'You must set different mo_class with qp set_mo_class'
+    print*, 'If you want to kick all the orbital:'
+    print*, 'qp set Orbital_optimization security_mo_class false'
+    print*, ''
+    print*, 'abort'
+
+    call abort
+  
+  endif
+  
+  call apply_pre_rotation
+  
+end
+#+END_SRC