Force MOs to be on axes. Nice for atoms

bin/qp_plugins

@@ -97,7 +97,7 @@ end
 
 def get_repositories():
     l_result = [f for f in os.listdir(QP_PLUGINS) \
-                  if f not in [".gitignore", "local"] ]
+                  if f not in [".gitignore", "local", "README.rst"] ]
     return sorted(l_result)
 
 

src/ao_one_e_ints/ao_ortho_canonical.irp.f

@@ -138,6 +138,8 @@ END_PROVIDER
     deallocate(S)
   endif
 
+  FREE ao_overlap
+
 END_PROVIDER
 
 BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)]

src/mo_basis/utils.irp.f

@@ -228,7 +228,11 @@ subroutine mo_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,eig,label)
   call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef,size(mo_coef,1))
 
   do i=1,m
-    eig(i) = D(i)
+    if (eig(i) > 1.d-20) then
+      eig(i) = D(i)
+    else
+      eig(i) = 0.d0
+    endif
   enddo
 
   deallocate(A,mo_coef_new,U,Vt,D)

src/scf_utils/diagonalize_fock.irp.f

@@ -47,7 +47,7 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num)
        do j = 1, n_core_orb
          jorb = list_core(j)
          F(iorb,jorb) = 0.d0
-         F(jorb,iorb) = 0.d0                                                                                                                 
+         F(jorb,iorb) = 0.d0
        enddo
      enddo
    endif

src/scf_utils/roothaan_hall_scf.irp.f

@@ -13,9 +13,9 @@ END_DOC
   integer                        :: iteration_SCF,dim_DIIS,index_dim_DIIS
 
   logical                        :: converged
-  integer                        :: i,j
+  integer                        :: i,j,m
   logical, external              :: qp_stop
-  double precision, allocatable :: mo_coef_save(:,:)
+  double precision, allocatable :: mo_coef_save(:,:), S(:,:)
 
   PROVIDE ao_md5 mo_occ level_shift
 
@@ -208,9 +208,29 @@ END_DOC
       size(Fock_matrix_mo,2),mo_label,1,.true.)
    call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10)
    call orthonormalize_mos
-   call save_mos
   endif
 
+
+  ! Identify degenerate MOs and force them on the axes
+  allocate(S(ao_num,ao_num))
+  i=1
+  do while (i<mo_num)
+    j=i
+    m=1
+    do while ( (j<mo_num).and.(fock_matrix_diag_mo(j+1)-fock_matrix_diag_mo(i) < 1.d-8) )
+      j += 1
+      m += 1
+    enddo
+    if (m>1) then
+      call dgemm('N','T',ao_num,ao_num,m,1.d0,mo_coef(1,i),size(mo_coef,1),mo_coef(1,i),size(mo_coef,1),0.d0,S,size(S,1))
+      call pivoted_cholesky( S, m, -1.d0, ao_num, mo_coef(1,i))
+    endif
+    i = j+1
+  enddo
+
+
+  call save_mos
+
   call write_double(6, Energy_SCF, 'SCF energy')
 
   call write_time(6)

src/trexio/export_trexio_routines.irp.f

@@ -59,7 +59,59 @@ subroutine export_trexio(update,full_path)
   enddo
   call ezfio_set_trexio_trexio_file(trexio_filename)
 
-    
+
+
+! ------------------------------------------------------------------------------
+
+! Metadata
+! --------
+
+  integer :: code_num, author_num
+  character*(64) :: code(100), author(100), user
+  character*(64), parameter :: qp2_code = "QuantumPackage"
+
+  call getenv("USER",user)
+  do k=1,N_states
+    rc = trexio_read_metadata_code_num(f(k), code_num)
+    if (rc == TREXIO_ATTR_MISSING) then
+      i = 1
+      code(:) = ""
+    else
+      rc = trexio_read_metadata_code(f(k), code, 64)
+      do i=1, code_num
+        if (trim(code(i)) == trim(qp2_code)) then
+          exit
+        endif
+      enddo
+    endif
+    if (i == code_num+1) then
+      code(i) = qp2_code
+      rc = trexio_write_metadata_code_num(f(k), i)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+      rc = trexio_write_metadata_code(f(k), code, 64)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    endif
+
+    rc = trexio_read_metadata_author_num(f(k), author_num)
+    if (rc == TREXIO_ATTR_MISSING) then
+      i = 1
+      author(:) = ""
+    else
+      rc = trexio_read_metadata_author(f(k), author, 64)
+      do i=1, author_num
+        if (trim(author(i)) == trim(user)) then
+          exit
+        endif
+      enddo
+    endif
+    if (i == author_num+1) then
+      author(i) = user
+      rc = trexio_write_metadata_author_num(f(k), i)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+      rc = trexio_write_metadata_author(f(k), author, 64)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    endif
+  enddo
 
 ! ------------------------------------------------------------------------------