Merge pull request #9 from pfloos/utest

Utest

.gitignore

@@ -1,2 +1,5 @@
 *.o
 *.
+__pycache__
+
+.ninja_deps

src/GT/GTpp_excitation_density.f90

@@ -135,131 +135,132 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
     dim_1 = (nBas - nO) * (nBas - nO - 1) / 2
     dim_2 = nO * (nO - 1) / 2
 
-    allocate(ERI_1(nBas,nBas,dim_1), ERI_2(nBas,nBas,dim_2))
+    if((dim_1 .eq. 0) .or. (dim_2 .eq. 0)) then
-    ERI_1 = 0.d0
-    ERI_2 = 0.d0
 
-    !$OMP PARALLEL DEFAULT(NONE)                     &
+      !$OMP PARALLEL DEFAULT(NONE)                                         &
-    !$OMP          PRIVATE(p, q, c, d, cd, k, l, kl) &
+      !$OMP          PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
-    !$OMP          SHARED(nC, nBas, nR, nO, ERI_1, ERI_2, ERI)
+      !$OMP          SHARED(nC, nBas, nR, nO, rho1, rho2, ERI, X1, Y1, X2, Y2)
-    !$OMP DO COLLAPSE(2)
+      !$OMP DO COLLAPSE(2)
-    do q = nC+1, nBas-nR
+      do q = nC+1, nBas-nR
-      do p = nC+1, nBas-nR
+        do p = nC+1, nBas-nR
-        cd = 0
+   
-        do c = nO+1, nBas-nR
+          ab = 0
-          do d = c+1, nBas-nR
+  
-            cd = cd + 1
+          do a = nO+1, nBas-nR
-            ERI_1(p,q,cd) = ERI(p,q,c,d) - ERI(p,q,d,c)
+            do b = a+1, nBas-nR
-          enddo
+  
-        enddo
+            ab = ab + 1
-        kl = 0
+   
-        do k = nC+1, nO
+            cd = 0
-          do l = k+1, nO
+            do c = nO+1, nBas-nR
-            kl = kl + 1
+              do d = c+1, nBas-nR
-            ERI_2(p,q,kl) = ERI(p,q,k,l) - ERI(p,q,l,k)
+  
-          end do
+                cd = cd + 1
-        end do
+  
-      enddo
+                rho1(p,q,ab) = rho1(p,q,ab) & 
-    enddo
+                             + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) 
+              end do ! d
+            end do ! c
+   
+            kl = 0
+            do k = nC+1, nO
+              do l = k+1, nO
+  
+                kl = kl + 1
+  
+                rho1(p,q,ab) = rho1(p,q,ab) & 
+                             + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) 
+              end do ! l
+            end do ! k 
+          end do ! b
+        end do ! a
+   
+        ij = 0
+        do i = nC+1, nO
+          do j = i+1, nO
+  
+            ij = ij + 1
+            
+            cd = 0
+            
+            do c = nO+1, nBas-nR
+              do d = c+1, nBas-nR
+            
+                cd = cd + 1
+            
+                rho2(p,q,ij) = rho2(p,q,ij) & 
+                             + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) 
+              end do ! d
+            end do ! c
+            
+            kl = 0
+            do k = nC+1, nO
+              do l = k+1, nO
+            
+                kl = kl + 1
+            
+                rho2(p,q,ij) = rho2(p,q,ij) & 
+                             + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) 
+              end do ! l
+            end do ! k
+          end do ! j
+        end do ! i
+      end do ! p
+    end do ! q
     !$OMP END DO
     !$OMP END PARALLEL
 
-    call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, &
+    else
-               ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, &
-               0.d0, rho1(1,1,1), nBas*nBas)
 
-    call dgemm("N", "N", nBas*nBas, dim_1, dim_2, 1.d0, &
+      allocate(ERI_1(nBas,nBas,dim_1), ERI_2(nBas,nBas,dim_2))
-               ERI_2(1,1,1), nBas*nBas, Y1(1,1), dim_2, &
+      ERI_1 = 0.d0
-               1.d0, rho1(1,1,1), nBas*nBas)
+      ERI_2 = 0.d0
   
-    call dgemm("N", "N", nBas*nBas, dim_2, dim_1, 1.d0, &
+      !$OMP PARALLEL DEFAULT(NONE)                     &
-               ERI_1(1,1,1), nBas*nBas, X2(1,1), dim_1, &
+      !$OMP          PRIVATE(p, q, c, d, cd, k, l, kl) &
-               0.d0, rho2(1,1,1), nBas*nBas)
+      !$OMP          SHARED(nC, nBas, nR, nO, ERI_1, ERI_2, ERI)
+      !$OMP DO COLLAPSE(2)
+      do q = nC+1, nBas-nR
+        do p = nC+1, nBas-nR
+          cd = 0
+          do c = nO+1, nBas-nR
+            do d = c+1, nBas-nR
+              cd = cd + 1
+              ERI_1(p,q,cd) = ERI(p,q,c,d) - ERI(p,q,d,c)
+            enddo
+          enddo
+          kl = 0
+          do k = nC+1, nO
+            do l = k+1, nO
+              kl = kl + 1
+              ERI_2(p,q,kl) = ERI(p,q,k,l) - ERI(p,q,l,k)
+            end do
+          end do
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
   
-    call dgemm("N", "N", nBas*nBas, dim_2, dim_2, 1.d0, &
+      call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, &
-               ERI_2(1,1,1), nBas*nBas, Y2(1,1), dim_2, &
+                 ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, &
-               1.d0, rho2(1,1,1), nBas*nBas)
+                 0.d0, rho1(1,1,1), nBas*nBas)
   
-    deallocate(ERI_1, ERI_2)
+      call dgemm("N", "N", nBas*nBas, dim_1, dim_2, 1.d0, &
+                 ERI_2(1,1,1), nBas*nBas, Y1(1,1), dim_2, &
+                 1.d0, rho1(1,1,1), nBas*nBas)
   
+      call dgemm("N", "N", nBas*nBas, dim_2, dim_1, 1.d0, &
+                 ERI_1(1,1,1), nBas*nBas, X2(1,1), dim_1, &
+                 0.d0, rho2(1,1,1), nBas*nBas)
   
-!    !$OMP PARALLEL DEFAULT(NONE)                                         &
+      call dgemm("N", "N", nBas*nBas, dim_2, dim_2, 1.d0, &
-!    !$OMP          PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
+                 ERI_2(1,1,1), nBas*nBas, Y2(1,1), dim_2, &
-!    !$OMP          SHARED(nC, nBas, nR, nO, rho1, rho2, ERI, X1, Y1, X2, Y2)
+                 1.d0, rho2(1,1,1), nBas*nBas)
-!    !$OMP DO COLLAPSE(2)
-!    do q = nC+1, nBas-nR
-!      do p = nC+1, nBas-nR
-! 
-!        ab = 0
-!
-!        do a = nO+1, nBas-nR
-!          do b = a+1, nBas-nR
-!
-!          ab = ab + 1
-! 
-!          cd = 0
-!          do c = nO+1, nBas-nR
-!            do d = c+1, nBas-nR
-!
-!              cd = cd + 1
-!
-!              rho1(p,q,ab) = rho1(p,q,ab) & 
-!                           + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) 
-!            end do ! d
-!          end do ! c
-! 
-!          kl = 0
-!          do k = nC+1, nO
-!            do l = k+1, nO
-!
-!              kl = kl + 1
-!
-!              rho1(p,q,ab) = rho1(p,q,ab) & 
-!                           + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) 
-!            end do ! l
-!          end do ! k
-! 
-!        end do ! b
-!      end do ! a
-! 
-!      ij = 0
-!      do i = nC+1, nO
-!        do j = i+1, nO
-!
-!          ij = ij + 1
-!          
-!          cd = 0
-!          
-!          do c = nO+1, nBas-nR
-!            do d = c+1, nBas-nR
-!          
-!              cd = cd + 1
-!          
-!              rho2(p,q,ij) = rho2(p,q,ij) & 
-!                           + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) 
-!            end do ! d
-!          end do ! c
-!          
-!          kl = 0
-!          do k = nC+1, nO
-!            do l = k+1, nO
-!          
-!              kl = kl + 1
-!          
-!              rho2(p,q,ij) = rho2(p,q,ij) & 
-!                           + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) 
-!            end do ! l
-!          end do ! k
-!
-!        end do ! j
-!      end do ! i
-!
-!    end do ! p
-!  end do ! q
-!  !$OMP END DO
-!  !$OMP END PARALLEL
   
-  end if
+      deallocate(ERI_1, ERI_2)
+
+    endif
+  endif
 
 !----------------------------------------------
 ! alpha-beta block
@@ -270,125 +271,126 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
     dim_1 = (nBas - nO) * (nBas - nO)
     dim_2 = nO * nO
 
-    allocate(ERI_1(nBas,nBas,dim_1), ERI_2(nBas,nBas,dim_2))
+    if((dim_1 .eq. 0) .or. (dim_2 .eq. 0)) then
-    ERI_1 = 0.d0
-    ERI_2 = 0.d0
 
-    !$OMP PARALLEL DEFAULT(NONE)                     &
+      !$OMP PARALLEL DEFAULT(NONE)                                         &
-    !$OMP          PRIVATE(p, q, c, d, cd, k, l, kl) &
+      !$OMP          PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
-    !$OMP          SHARED(nC, nBas, nR, nO, ERI_1, ERI_2, ERI)
+      !$OMP          SHARED(nC, nBas, nR, nO, rho1, rho2, ERI, X1, Y1, X2, Y2)
-    !$OMP DO COLLAPSE(2)
+      !$OMP DO COLLAPSE(2)
-    do q = nC+1, nBas-nR
+       
-      do p = nC+1, nBas-nR
+      do q = nC+1, nBas-nR
-        cd = 0
+        do p = nC+1, nBas-nR
-        do c = nO+1, nBas-nR
+   
-          do d = nO+1, nBas-nR
+          ab = 0
-            cd = cd + 1
+          do a = nO+1, nBas-nR
-            ERI_1(p,q,cd) = ERI(p,q,c,d)
+            do b = nO+1, nBas-nR
-          enddo
+  
-        enddo
+              ab = ab + 1
-        kl = 0
+   
-        do k = nC+1, nO
+              cd = 0
-          do l = nC+1, nO
+              do c = nO+1, nBas-nR
-            kl = kl + 1
+                do d = nO+1, nBas-nR
-            ERI_2(p,q,kl) = ERI(p,q,k,l)
+  
+                  cd = cd + 1
+  
+                  rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,c,d)*X1(cd,ab) 
+                end do
+              end do
+   
+              kl = 0
+              do k = nC+1, nO
+                do l = nC+1, nO
+  
+                  kl = kl + 1
+  
+                  rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,k,l)*Y1(kl,ab) 
+                end do
+              end do
+            end do
+          end do
+   
+          ij = 0
+          do i = nC+1, nO
+            do j = nC+1, nO
+  
+              ij = ij + 1
+               
+              cd = 0
+              do c = nO+1, nBas-nR
+                do d = nO+1, nBas-nR
+  
+                  cd = cd + 1
+  
+                  rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,c,d)*X2(cd,ij) 
+                end do
+              end do
+   
+              kl = 0
+              do k = nC+1, nO
+                do l = nC+1, nO
+  
+                  kl = kl + 1
+  
+                  rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,k,l)*Y2(kl,ij) 
+                end do
+              end do
+            end do
           end do
         end do
+      end do
+      !$OMP END DO
+      !$OMP END PARALLEL
+
+    else
+
+      allocate(ERI_1(nBas,nBas,dim_1), ERI_2(nBas,nBas,dim_2))
+      ERI_1 = 0.d0
+      ERI_2 = 0.d0
+  
+      !$OMP PARALLEL DEFAULT(NONE)                     &
+      !$OMP          PRIVATE(p, q, c, d, cd, k, l, kl) &
+      !$OMP          SHARED(nC, nBas, nR, nO, ERI_1, ERI_2, ERI)
+      !$OMP DO COLLAPSE(2)
+      do q = nC+1, nBas-nR
+        do p = nC+1, nBas-nR
+          cd = 0
+          do c = nO+1, nBas-nR
+            do d = nO+1, nBas-nR
+              cd = cd + 1
+              ERI_1(p,q,cd) = ERI(p,q,c,d)
+            enddo
+          enddo
+          kl = 0
+          do k = nC+1, nO
+            do l = nC+1, nO
+              kl = kl + 1
+              ERI_2(p,q,kl) = ERI(p,q,k,l)
+            end do
+          end do
+        enddo
       enddo
-    enddo
+      !$OMP END DO
-    !$OMP END DO
+      !$OMP END PARALLEL
-    !$OMP END PARALLEL
   
-    call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, &
+      call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, &
-               ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, &
+                 ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, &
-               0.d0, rho1(1,1,1), nBas*nBas)
+                 0.d0, rho1(1,1,1), nBas*nBas)
   
-    call dgemm("N", "N", nBas*nBas, dim_1, dim_2, 1.d0, &
+      call dgemm("N", "N", nBas*nBas, dim_1, dim_2, 1.d0, &
-               ERI_2(1,1,1), nBas*nBas, Y1(1,1), dim_2, &
+                 ERI_2(1,1,1), nBas*nBas, Y1(1,1), dim_2, &
-               1.d0, rho1(1,1,1), nBas*nBas)
+                 1.d0, rho1(1,1,1), nBas*nBas)
   
-    call dgemm("N", "N", nBas*nBas, dim_2, dim_1, 1.d0, &
+      call dgemm("N", "N", nBas*nBas, dim_2, dim_1, 1.d0, &
-               ERI_1(1,1,1), nBas*nBas, X2(1,1), dim_1, &
+                 ERI_1(1,1,1), nBas*nBas, X2(1,1), dim_1, &
-               0.d0, rho2(1,1,1), nBas*nBas)
+                 0.d0, rho2(1,1,1), nBas*nBas)
   
-    call dgemm("N", "N", nBas*nBas, dim_2, dim_2, 1.d0, &
+      call dgemm("N", "N", nBas*nBas, dim_2, dim_2, 1.d0, &
-               ERI_2(1,1,1), nBas*nBas, Y2(1,1), dim_2, &
+                 ERI_2(1,1,1), nBas*nBas, Y2(1,1), dim_2, &
-               1.d0, rho2(1,1,1), nBas*nBas)
+                 1.d0, rho2(1,1,1), nBas*nBas)
   
-    deallocate(ERI_1, ERI_2)
+      deallocate(ERI_1, ERI_2)
 
-
+    endif
-!    !$OMP PARALLEL DEFAULT(NONE)                                         &
+  endif
-!    !$OMP          PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
-!    !$OMP          SHARED(nC, nBas, nR, nO, rho1, rho2, ERI, X1, Y1, X2, Y2)
-!    !$OMP DO COLLAPSE(2)
-!     
-!    do q = nC+1, nBas-nR
-!      do p = nC+1, nBas-nR
-! 
-!        ab = 0
-!        do a = nO+1, nBas-nR
-!          do b = nO+1, nBas-nR
-!
-!            ab = ab + 1
-! 
-!            cd = 0
-!            do c = nO+1, nBas-nR
-!              do d = nO+1, nBas-nR
-!
-!                cd = cd + 1
-!
-!                rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,c,d)*X1(cd,ab) 
-!              end do
-!            end do
-! 
-!            kl = 0
-!            do k = nC+1, nO
-!              do l = nC+1, nO
-!
-!                kl = kl + 1
-!
-!                rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,k,l)*Y1(kl,ab) 
-!              end do
-!            end do
-! 
-!          end do
-!        end do
-! 
-!        ij = 0
-!        do i = nC+1, nO
-!          do j = nC+1, nO
-!
-!            ij = ij + 1
-!             
-!            cd = 0
-!            do c = nO+1, nBas-nR
-!              do d = nO+1, nBas-nR
-!
-!                cd = cd + 1
-!
-!                rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,c,d)*X2(cd,ij) 
-!              end do
-!            end do
-! 
-!            kl = 0
-!            do k = nC+1, nO
-!              do l = nC+1, nO
-!
-!                kl = kl + 1
-!
-!                rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,k,l)*Y2(kl,ij) 
-!              end do
-!            end do
-! 
-!          end do
-!        end do
-!       
-!      end do
-!    end do
-!    !$OMP END DO
-!    !$OMP END PARALLEL
-
-  end if
 
 end subroutine 

src/GW/evRGW.f90

@@ -209,7 +209,9 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 ! Cumulant expansion !
 !--------------------!
 
-  call RGWC(dotest,nBas,nC,nO,nR,nS,Om,rho,eGW,Z)
+  ! TODO
+  !call RGWC(dotest, eta, nBas, nC, nO, nV, nR, nS, Om, rho, eHF, eGW, eGW, Z)
+  call RGWC(dotest, eta, nBas, nC, nO, nV, nR, nS, Om, rho, eHF, eHF, eGW, Z)
 
 ! Deallocate memory
 

src/LR/ppLR.f90

@@ -72,38 +72,43 @@ subroutine ppLR(TDA, nOO, nVV, Bpp, Cpp, Dpp, Om1, X1, Y1, Om2, X2, Y2, EcRPA)
     M(    1:nVV    ,nVV+1:nOO+nVV) = -           Bpp(1:nVV,1:nOO)
     M(nVV+1:nOO+nVV,    1:nVV)     = + transpose(Bpp(1:nVV,1:nOO))
 
-    !! Diagonalize the p-p matrix
+    if((nOO .eq. 0) .or. (nVV .eq. 0)) then
-    !if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z)
-    !! Split the various quantities in p-p and h-h parts
-    !call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2)
 
+      ! Diagonalize the p-p matrix
+      if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z)
+      ! Split the various quantities in p-p and h-h parts
+      call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2)
 
-    thr_d   = 1d-6   ! to determine if     diagonal elements of L.T x R are close enouph to 1
+    else
-    thr_nd  = 1d-6   ! to determine if non-diagonal elements of L.T x R are close enouph to 1
-    thr_deg = 1d-8   ! to determine if two eigenvectors are degenerate or not
-    imp_bio = .True. ! impose bi-orthogonality
-    verbose = .False.
-    call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
 
-    do i = 1, nOO
+        thr_d   = 1d-6   ! to determine if     diagonal elements of L.T x R are close enouph to 1
-      Om2(i) = Om(i)
+        thr_nd  = 1d-6   ! to determine if non-diagonal elements of L.T x R are close enouph to 1
-      do j = 1, nVV
+        thr_deg = 1d-8   ! to determine if two eigenvectors are degenerate or not
-        X2(j,i) = Z(j,i)
+        imp_bio = .True. ! impose bi-orthogonality
-      enddo
+        verbose = .False.
-      do j = 1, nOO
+        call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
-        Y2(j,i) = Z(nVV+j,i)
-      enddo
-    enddo
     
-    do i = 1, nVV
+        do i = 1, nOO
-      Om1(i) = Om(nOO+i)
+          Om2(i) = Om(i)
-      do j = 1, nVV
+          do j = 1, nVV
-        X1(j,i) = M(j,nOO+i)
+            X2(j,i) = Z(j,i)
-      enddo
+          enddo
-      do j = 1, nOO
+          do j = 1, nOO
-        Y1(j,i) = M(nVV+j,nOO+i)
+            Y2(j,i) = Z(nVV+j,i)
-      enddo
+          enddo
-    enddo
+        enddo
+    
+        do i = 1, nVV
+          Om1(i) = Om(nOO+i)
+          do j = 1, nVV
+            X1(j,i) = M(j,nOO+i)
+          enddo
+          do j = 1, nOO
+            Y1(j,i) = M(nVV+j,nOO+i)
+          enddo
+        enddo
+
+    endif
 
   end if
 

src/test/check_test_value.f90

@@ -9,12 +9,12 @@ subroutine check_test_value(branch)
 ! Local variables
 
   character(len=30)                  :: description
-  double precision                   :: value
+  double precision                   :: val
   double precision                   :: reference
   character(len=15)                  :: answer
 
   logical                            :: failed
-  double precision,parameter         :: cutoff = 1d-10
+  double precision,parameter         :: thresh = 1d-10
 
 ! Output variables
 
@@ -45,19 +45,19 @@ subroutine check_test_value(branch)
   do
 
     read(11,'(A30)',end=11) description
-    read(11,'(F20.15)',end=11) value
+    read(11,'(F20.15)',end=11) val
 
     read(12,*,end=12) 
     read(12,'(F20.15)',end=12) reference
 
-    if(abs(value-reference) < cutoff) then
+    if(dabs(val-reference)/(1d-15+dabs(reference)) < thresh) then
       answer = '.......... :-)'
     else
       answer = '.......... :-( '
       failed = .true.
     end if
     write(*,'(1X,A1,1X,A30,1X,A1,1X,3F15.10,1X,A1,1X,A15,1X,A1)') & 
-      '|',description,'|',value,reference,abs(value-reference),'|',answer,'|'
+      '|',description,'|',val,reference,abs(val-reference),'|',answer,'|'
 
   end do
 

src/test/dump_test_value.f90

@@ -1,4 +1,4 @@
-subroutine dump_test_value(branch,description,value)
+subroutine dump_test_value(branch, description, val)
 
   implicit none
 
@@ -7,7 +7,7 @@ subroutine dump_test_value(branch,description,value)
   character(len=1),intent(in)        :: branch
   character(len=*),intent(in)        :: description
 
-  double precision,intent(in)        :: value
+  double precision,intent(in)        :: val
  
 ! Local variables
 
@@ -15,18 +15,19 @@ subroutine dump_test_value(branch,description,value)
 
   if(branch == 'R') then 
 
-    write(11,*) trim(description)
+    !write(1231597, '(A, ": ", F20.15)') '"' // trim(description) // '"', val
-    write(11,'(F20.15)') value
+    write(1231597, *) trim(description)
+    write(1231597, '(F20.15)') val
 
   elseif(branch == 'U') then 
 
-    write(12,*) trim(description)
+    write(1232584,*) trim(description)
-    write(12,'(F20.15)') value
+    write(1232584,'(F20.15)') val
 
   elseif(branch == 'G') then
 
-    write(13,*) trim(description)
+    write(1234181,*) trim(description)
-    write(13,'(F20.15)') value
+    write(1234181,'(F20.15)') val
 
   else
 

src/test/init_test.f90

@@ -12,10 +12,10 @@ subroutine init_test(doRtest,doUtest,doGtest)
 
 ! Output variables
 
-  if(doRtest) open(unit=11,file='test/Rtest.dat')
+  if(doRtest) open(unit=1231597, file='test/Rtest.dat')
 
-  if(doUtest) open(unit=12,file='test/Utest.dat')
+  if(doUtest) open(unit=1232584, file='test/Utest.dat')
 
-  if(doGtest) open(unit=13,file='test/Gtest.dat')
+  if(doGtest) open(unit=1234181, file='test/Gtest.dat')
 
 end subroutine

src/test/stop_test.f90

@@ -12,10 +12,10 @@ subroutine stop_test(doRtest,doUtest,doGtest)
 
 ! Output variables
 
-  if(doRtest) close(unit=11)
+  if(doRtest) close(unit=1231597)
 
-  if(doUtest) close(unit=12)
+  if(doUtest) close(unit=1231597)
 
-  if(doGtest) close(unit=13)
+  if(doGtest) close(unit=1234181)
 
 end subroutine

src/utils/non_sym_diag.f90

@@ -3,7 +3,7 @@
 
 subroutine diagonalize_nonsym_matrix(N, A, L, e_re, thr_d, thr_nd, thr_deg, imp_bio, verbose)
 
-  ! Diagonalize a non-symmetric matrix
+  ! Diagonalize a non-symmetric matrix A
   ! 
   ! Output
   !   right-eigenvectors are saved in A
@@ -278,7 +278,7 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot)
 
     do i = 1, m
       if(S(i,i) <= 0.d0) then
-        print *, ' overap negative'
+        print *, ' negative overlap !'
         print *, i, S(i,i)
         exit
       endif

test/export_tobench.py

@@ -0,0 +1,26 @@
+
+import sys
+
+
+def read_quantities_from_file(filename):
+    quantities = {}
+
+    with open(filename, 'r') as file:
+        lines = file.readlines()
+        for i in range(0, len(lines), 2):
+            # Remove any leading or trailing whitespace/newline characters
+            quantity_name = lines[i].strip()
+            quantity_value = float(lines[i+1].strip())
+            quantities[quantity_name] = quantity_value
+
+    return quantities
+
+def print_quantities(quantities):
+    for key, value in quantities.items():
+        print(f'"{key}": {value},')
+
+filename = sys.argv[1]
+
+quantities = read_quantities_from_file(filename)
+print_quantities(quantities)
+

tests/.gitignore

@@ -0,0 +1,7 @@
+
+FeatherBench.db
+FeatherBench.json
+
+*.xyz
+work
+

tests/create_database.py

@@ -0,0 +1,52 @@
+
+import argparse
+
+from molecule import save_molecules_to_json, load_molecules_from_json
+from molecule import create_database, add_molecule_to_db, remove_database
+
+from feather_bench import FeatherBench
+
+
+parser = argparse.ArgumentParser(description="Benchmark Data Sets")
+
+parser.add_argument(
+    '-s', '--set_type',
+    choices=['light', 'medium', 'heavy'],
+    default='light',
+    help="Specify the type of data set: light (default), medium, or heavy."
+)
+
+
+args = parser.parse_args()
+
+if args.set_type == 'light':
+    bench = 'FeatherBench'
+    bench_title = "\n\nSelected Light Benchmark: {}\n\n".format(bench)
+elif args.set_type == 'medium':
+    bench = 'BalanceBench'
+    bench_title = "\n\nSelected Medium Benchmark: {}\n\n".format(bench)
+elif args.set_type == 'heavy':
+    bench = 'TitanBench'
+    bench_title = "\n\nSelected Heavy Benchmark: {}\n\n".format(bench)
+else:
+    bench_title = "\n\nSelected Light Benchmark: {}\n\n".format(bench)
+
+
+db_name = '{}.db'.format(bench)
+
+
+# Save molecules to JSON
+#save_molecules_to_json(FeatherBench, 'FeatherBench.json')
+
+# Load molecules from JSON
+#loaded_molecules = load_molecules_from_json('FeatherBench.json')
+#print(loaded_molecules)
+
+#remove_database(db_name)
+
+create_database(db_name)
+for molecule in FeatherBench:
+    add_molecule_to_db(db_name, molecule)
+
+
+

tests/feather_bench.py

@@ -0,0 +1,113 @@
+
+from molecule import Molecule
+
+
+He = Molecule(
+    name="He",
+    multiplicity=1,
+    geometry=[
+        {"element": "He", "x": 0.0, "y": 0.0, "z": 0.0}
+    ],
+    properties={
+        "properties_rhf":{
+            "6-31g": {
+                "RHF energy": -2.855160426154444,
+                "RHF HOMO energy": -0.914126628640145,
+                "RHF LUMO energy": 1.399859335255765,
+                "RHF dipole moment": 0.0,
+                "MP2 correlation energy": -0.011200122909934,
+                "CCD correlation energy": -0.014985063116,
+                "CCSD correlation energy": -0.015001711549092,
+                "drCCD correlation energy": -0.01884537385338,
+                "rCCD correlation energy": -0.016836322809386,
+                "crCCD correlation energy": 0.008524676641474,
+                "lCCD correlation energy": -0.00808242082105,
+                "CIS singlet excitation energy": 1.911193619991987,
+                "CIS triplet excitation energy": 1.455852629458543,
+                "phRPA correlation energy": -0.018845374128748,
+                "phRPAx correlation energy": -0.015760565120758,
+                "crRPA correlation energy": -0.008868581132249,
+                "ppRPA correlation energy": -0.008082420814972,
+                "G0F2 correlation energy": -0.011438430540104,
+                "G0F2 HOMO energy": -0.882696116274599,
+                "G0F2 LUMO energy": 1.383080391842522,
+                "G0W0 correlation energy": -0.019314094399372,
+                "G0W0 HOMO energy": -0.87053388021722,
+                "G0W0 LUMO energy": 1.377171287041735,
+                "evGW correlation energy": -0.019335511771337,
+                "evGW HOMO energy": -0.868460640984803,
+                "evGW LUMO energy": 1.376287581502582,
+                "G0T0pp correlation energy": -0.008161908540634,
+                "G0T0pp HOMO energy": -0.898869172597701,
+                "G0T0pp LUMO energy": 1.383928087417952,
+            }
+        },
+        "properties_uhf":{
+            "6-31g": {
+            }
+        },
+        "properties_ghf":{
+            "6-31g": {
+            }
+        },
+    }
+)
+
+# ---
+
+H2O = Molecule(
+    name="H2O",
+    multiplicity=1,
+    geometry=[
+        {"element": "O", "x":  0.0000, "y": 0.0000, "z": 0.0000},
+        {"element": "H", "x":  0.7571, "y": 0.0000, "z": 0.5861},
+        {"element": "H", "x": -0.7571, "y": 0.0000, "z": 0.5861}
+    ],
+    properties={
+        "properties_rhf":{
+            "cc-pvdz": {
+                "RHF energy": -85.21935817501823,
+                "RHF HOMO energy": -0.493132793449897,
+                "RHF LUMO energy": 0.185534869842355,
+                "RHF dipole moment": 0.233813698748474,
+                "MP2 correlation energy": -0.203978216774657,
+                "CCD correlation energy": -0.212571260121257,
+                "CCSD correlation energy": -0.213302190845899,
+                "drCCD correlation energy": -0.231281853419338,
+                "rCCD correlation energy": -0.277238348710547,
+                "crCCD correlation energy": 0.18014617422324,
+                "lCCD correlation energy": -0.15128653432796,
+                "CIS singlet excitation energy": 0.338828950934568,
+                "CIS triplet excitation energy": 0.304873339484139,
+                "phRPA correlation energy": -0.231281866582435,
+                "phRPAx correlation energy": -0.310796738307943,
+                "crRPA correlation energy": -0.246289801609294,
+                "ppRPA correlation energy": -0.151286536255888,
+                "G0F2 correlation energy": -0.217807591229668,
+                "G0F2 HOMO energy": -0.404541451101377,
+                "G0F2 LUMO energy": 0.16650398400197,
+                "G0W0 correlation energy": -0.23853664665404,
+                "G0W0 HOMO energy": -0.446828623007469,
+                "G0W0 LUMO energy": 0.173026609033024,
+                "evGW correlation energy": -0.239414217281308,
+                "evGW HOMO energy": -0.443076613314424,
+                "evGW LUMO energy": 0.172691758111392,
+                "G0T0pp correlation energy": -0.156214864467344,
+                "G0T0pp HOMO energy": -0.452117482732615,
+                "G0T0pp LUMO energy": 0.16679206983464,
+            }
+        }
+    }
+)
+
+# ---
+
+FeatherBench = [
+    He, 
+    H2O
+]
+
+
+
+
+

tests/inp/methods.RHF

@@ -0,0 +1,22 @@
+# RHF UHF GHF ROHF
+  T   F   F   F
+# MP2 MP3
+  T   T
+# CCD pCCD DCD CCSD CCSD(T)
+  T   F    F   T    F
+# drCCD rCCD crCCD lCCD
+  T     T    T     T
+# CIS CIS(D) CID CISD FCI
+  T   F      F   F    F
+# phRPA phRPAx crRPA ppRPA
+  T     T      T     T 
+# G0F2 evGF2 qsGF2 ufGF2 G0F3 evGF3
+  T    F     F     F    F    F
+# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
+  T    T    F    F        F      F
+# G0T0pp evGTpp qsGTpp ufG0T0pp
+  T      F      F      F
+# G0T0eh evGTeh qsGTeh
+  F      F      F
+# Rtest Utest Gtest
+  T     F     F

tests/inp/options.RHF

@@ -0,0 +1,18 @@
+# HF: maxSCF thresh  DIIS guess mix   shift stab search
+      10000  0.0000001  5 1     0.0   0.0   F    F
+# MP: reg
+      F
+# CC: maxSCF thresh     DIIS
+      64     0.0000001  5
+# spin: TDA singlet triplet
+        F   T       T
+# GF: maxSCF thresh  DIIS lin eta renorm reg
+      256    0.00001 5    F   0.0 0      F
+# GW: maxSCF thresh  DIIS lin eta TDA_W reg
+      256   0.00001  5    F   0.0 F     F 
+# GT: maxSCF thresh  DIIS lin eta TDA_T reg
+      256    0.00001 5    F   0.0 F     F  
+# ACFDT: AC Kx  XBS
+         F  F   T
+# BSE: phBSE phBSE2 ppBSE dBSE dTDA
+       F     F      F     F    T     

tests/lunch_bench.py

@@ -0,0 +1,246 @@
+
+import time
+import threading
+import sys
+import os
+import shutil
+from pathlib import Path
+import subprocess
+import platform
+from datetime import datetime
+import argparse
+
+from molecule import get_molecules_from_db
+from molecule import generate_xyz
+from utils import print_col, stdout_col
+
+
+current_date = datetime.now()
+
+quack_root = os.getenv('QUACK_ROOT')
+
+# User Name
+user_name = os.getlogin()
+
+# Operating System
+os_name = platform.system()
+os_release = platform.release()
+os_version = platform.version()
+
+# CPU Information
+machine = platform.machine()
+processor = platform.processor()
+
+# System Architecture
+architecture = platform.architecture()[0]
+
+# Python Version
+python_version_full = platform.python_version_tuple()
+PYTHON_VERSION = "{}.{}".format(python_version_full[0], python_version_full[1])
+
+
+
+print(f"The current date and time is {current_date.strftime('%Y-%m-%d %H:%M:%S')}")
+print(f"User Name: {user_name}")
+print(f"Operating System: {os_name} {os_release} ({os_version})")
+print(f"CPU: {processor} ({machine})")
+print(f"System Architecture: {architecture}")
+print(f"QUACK_ROOT: {quack_root}")
+print(f"Python version: {python_version_full}\n\n")
+
+
+parser = argparse.ArgumentParser(description="Benchmark Data Sets")
+
+parser.add_argument(
+    '-s', '--set_type',
+    choices=['light', 'medium', 'heavy'],
+    default='light',
+    help="Specify the type of data set: light (default), medium, or heavy."
+)
+parser.add_argument(
+    '-t', '--thresh',
+    type=float,
+    default=1e-7,
+    help='Threshold for acceptable difference (default: 1e-8)'
+)
+    
+
+
+
+
+args = parser.parse_args()
+
+THRESH = args.thresh
+
+if args.set_type == 'light':
+    bench = 'FeatherBench'
+    bench_title = "\n\nSelected Light Benchmark: {}\n\n".format(bench)
+elif args.set_type == 'medium':
+    bench = 'BalanceBench'
+    bench_title = "\n\nSelected Medium Benchmark: {}\n\n".format(bench)
+elif args.set_type == 'heavy':
+    bench = 'TitanBench'
+    bench_title = "\n\nSelected Heavy Benchmark: {}\n\n".format(bench)
+else:
+    bench_title = "\n\nSelected Light Benchmark: {}\n\n".format(bench)
+
+print(bench_title.center(150, '-'))
+print("\n\n")
+
+# ---
+
+class Quack_Job:
+
+    def __init__(self, mol, multip, basis, geom, methd):
+        self.mol = mol
+        self.multip = multip
+        self.basis = basis
+        self.geom = geom
+        self.methd = methd
+
+    def prep_inp(self):
+
+        # geometry
+        generate_xyz(self.geom, filename="{}/mol/{}.xyz".format(quack_root, self.mol))
+
+        # input files
+        for inp in ["methods", "options"]:
+            inp_file = "{}.{}".format(inp, self.methd.upper())
+            if os.path.exists("inp/{}".format(inp_file)):
+                shutil.copy("{}/tests/inp/{}".format(quack_root, inp_file), 
+                            "{}/input/{}".format(quack_root, inp))
+            else:
+                print_col("File 'inp/{}' does not exist.".format(inp_file), "red")
+                sys.exit(1)
+
+    def run(self, work_path):
+
+        def display_spinner():
+            spinner = ['|', '/', '-', '\\']
+            idx = 0
+            while not done_event.is_set():
+                stdout_col(f'\r    Testing {self.methd} ({self.basis}) {spinner[idx]}', "cyan")
+                sys.stdout.flush()
+                idx = (idx + 1) % len(spinner)
+                time.sleep(0.05)
+            stdout_col(f'\r    Testing {self.methd} ({self.basis})    \n\n', "cyan")
+
+        done_event = threading.Event()
+        spinner_thread = threading.Thread(target=display_spinner)
+        spinner_thread.start()
+
+        try:
+    
+            os.chdir('..')
+            #print_col(f"      Starting QuAck..", "magenta")
+            #print_col(f"      $ cd ..", "magenta")
+    
+            command = [
+                'python{}'.format(PYTHON_VERSION), 'PyDuck.py',
+                '-x', '{}'.format(self.mol), 
+                '-b', '{}'.format(self.basis),
+                '-m', '{}'.format(self.multip)
+            ]
+            #print_col(f"      $ {' '.join(command)}", "magenta")
+    
+            file_out = "{}/{}/{}_{}_{}.out".format(work_path, self.methd, self.mol, self.multip, self.basis)
+            with open(file_out, 'w') as fobj:
+                result = subprocess.run(command, stdout=fobj, stderr=subprocess.PIPE, text=True)
+            if result.stderr:
+                print("Error output:", result.stderr)
+    
+            os.chdir('tests')
+            #print_col(f"      $ cd tests", "magenta")
+    
+        except Exception as e:
+
+            print_col(f"An error occurred: {str(e)}", "red")
+
+        finally:
+
+            done_event.set()
+            spinner_thread.join()
+
+    def check_data(self, data_ref):
+        filepath = '../test/Rtest.dat'
+        data_new = {}
+        try:
+            # read data_new
+            with open(filepath, 'r') as f:
+                lines = f.readlines()
+                for i in range(0, len(lines) - 1, 2):
+                    key = lines[i].strip()
+                    value = lines[i + 1].strip()
+                    data_new[key] = float(value)  # Convert value to float
+            
+            # Compare with data_ref
+            for key in data_ref:
+                if key not in data_new:
+                    print_col(f"        😐 {key} missing ⚠️ ", "yellow")
+                else:
+                    diff = abs(data_new[key] - data_ref[key]) / (1e-15 + abs(data_ref[key]))
+                    if(diff <= THRESH):
+                        print_col(f"        🙂 {key}", "green")
+                    else:
+                        print_col(f"        ☹️  {key}: ❌ {data_ref[key]} ≠ {data_new[key]}", "red")
+        except FileNotFoundError:
+            print_col(f"Error: The file '{filepath}' does not exist.", "red")
+            sys.exit(1)
+        except Exception as e:
+            print_col(f"An error occurred: {str(e)}", "red")
+            sys.exit(1)
+    
+
+# ---
+
+
+def main():
+
+    work_path = Path('{}/tests/work'.format(quack_root))
+    if not work_path.exists():
+        work_path.mkdir(parents=True, exist_ok=True)
+        print(f"Directory '{work_path}' created.\n")
+
+    for mol in molecules:
+
+        mol_name = mol.name
+        mol_mult = mol.multiplicity
+        mol_geom = mol.geometry
+        mol_data = mol.properties
+
+        print_col("  Molecule: {} (2S+1 = {})".format(mol_name, mol_mult), "blue")
+
+        for mol_prop_name, mol_prop_data in mol_data.items():
+
+            methd = mol_prop_name[len('properties_'):]
+
+            if(len(mol_prop_data) == 0):
+                continue
+
+            for basis_name, basis_data in mol_prop_data.items():
+
+                if(len(basis_data) == 0):
+                    continue
+
+                work_methd = Path('{}/{}'.format(work_path, methd))
+                if not work_methd.exists():
+                    work_methd.mkdir(parents=True, exist_ok=True)
+        
+                New_Quack_Job = Quack_Job(mol_name, mol_mult, basis_name, mol_geom, methd)
+                New_Quack_Job.prep_inp()
+                New_Quack_Job.run(work_path)
+                New_Quack_Job.check_data(basis_data)
+
+                print()
+            print()
+        print()
+
+    quit()
+
+        
+db_name = '{}.db'.format(bench)
+
+molecules = get_molecules_from_db(db_name)
+
+main()
+

tests/molecule.py

@@ -0,0 +1,150 @@
+
+import os
+import json
+import sqlite3
+
+from utils import print_col
+
+
+class Molecule:
+    def __init__(self, name, multiplicity, geometry, properties):
+        self.name = name
+        self.multiplicity = multiplicity
+        self.geometry = geometry
+        self.properties = properties
+
+    def to_dict(self):
+        return {
+            "name": self.name,
+            "multiplicity": self.multiplicity,
+            "geometry": self.geometry,
+            "properties": self.properties,
+        }
+
+    @staticmethod
+    def from_dict(data):
+        return Molecule(
+            name=data["name"],
+            multiplicity=data["multiplicity"],
+            geometry=data["geometry"],
+            properties=data["properties"]
+        )
+
+def save_molecules_to_json(molecules, filename):
+    with open(filename, 'w') as f:
+        json_data = [molecule.to_dict() for molecule in molecules]
+        json.dump(json_data, f, indent=4)
+
+def load_molecules_from_json(filename):
+    with open(filename, 'r') as f:
+        json_data = json.load(f)
+        return [Molecule.from_dict(data) for data in json_data]
+
+
+def create_database(db_name):
+    if os.path.exists(db_name):
+        conn = sqlite3.connect(db_name)
+        cursor = conn.cursor()
+        # Check if the table already exists
+        cursor.execute("SELECT name FROM sqlite_master WHERE type='table' AND name='molecules';")
+        table_exists = cursor.fetchone()
+        
+        if table_exists:
+            print_col(f"Database '{db_name}' already exists and table 'molecules' is already created.", "yellow")
+        else:
+            # Create the table if it does not exist
+            cursor.execute('''CREATE TABLE molecules
+                              (name TEXT, multiplicity INTEGER, geometry TEXT, properties TEXT)''')
+            conn.commit()
+            print_col(f"Table 'molecules' created in existing database '{db_name}' successfully.", "green")
+        conn.close()
+    else:
+        # Create the database and table
+        conn = sqlite3.connect(db_name)
+        cursor = conn.cursor()
+        cursor.execute('''CREATE TABLE molecules
+                          (name TEXT, multiplicity INTEGER, geometry TEXT, properties TEXT)''')
+        conn.commit()
+        conn.close()
+        print_col(f"Database '{db_name}' created and table 'molecules' added successfully.", "green")
+
+def add_molecule_to_db(db_name, molecule):
+
+    conn = sqlite3.connect(db_name)
+    cursor = conn.cursor()
+
+    # Convert geometry and properties to JSON strings
+    geometry_str = json.dumps(molecule.geometry)
+    energies_str = json.dumps(molecule.properties)
+
+    # Check if the molecule already exists
+    cursor.execute("SELECT COUNT(*) FROM molecules WHERE name = ?", (molecule.name,))
+    count = cursor.fetchone()[0]
+    
+    if count > 0:
+        print_col(f"Molecule '{molecule.name}' already exists in {db_name}.", "yellow")
+    else:
+        # Insert the molecule if it does not exist
+        cursor.execute("INSERT INTO molecules (name, multiplicity, geometry, properties) VALUES (?, ?, ?, ?)",
+                       (molecule.name, molecule.multiplicity, geometry_str, energies_str))
+        conn.commit()
+        print_col(f"'{molecule.name}' added to {db_name} successfully.", "green")
+    
+    conn.close()
+
+
+def remove_database(db_name):
+    if os.path.exists(db_name):
+        os.remove(db_name)
+        print_col(f"Database '{db_name}' removed successfully.", "red")
+    else:
+        print_col(f"Database '{db_name}' does not exist.", "red")
+
+def get_molecules_from_db(db_name):
+    conn = sqlite3.connect(db_name)
+    cursor = conn.cursor()
+    cursor.execute("SELECT name, multiplicity, geometry, properties FROM molecules")
+    rows = cursor.fetchall()
+    molecules = []
+    for row in rows:
+        name, multiplicity, geometry_str, energies_str = row
+        geometry = json.loads(geometry_str)
+        properties = json.loads(energies_str)
+        molecules.append(Molecule(name, multiplicity, geometry, properties))
+    conn.close()
+    return molecules
+
+
+def generate_xyz(elements, filename="output.xyz", verbose=False):
+    """
+    Generate an XYZ file from a list of elements.
+
+    Parameters:
+    elements (list): A list of dictionaries, where each dictionary represents
+                     an atom with its element and x, y, z coordinates.
+    filename (str): The name of the output XYZ file. Default is 'output.xyz'.
+    """
+
+    # Get the number of atoms
+    num_atoms = len(elements)
+
+    # Open the file in write mode
+    with open(filename, 'w') as f:
+        # Write the number of atoms
+        f.write(f"{num_atoms}\n")
+
+        # Write a comment line (can be left blank or customized)
+        f.write("XYZ file generated by generate_xyz function\n")
+
+        # Write the element and coordinates
+        for atom in elements:
+            element = atom['element']
+            x = atom['x']
+            y = atom['y']
+            z = atom['z']
+            f.write(f"{element} {x:.6f} {y:.6f} {z:.6f}\n")
+
+    if(verbose):
+        print(f"XYZ file '{filename}' generated successfully!")
+
+

tests/utils.py

@@ -0,0 +1,88 @@
+
+import sys
+
+
+
+
+def print_col(text, color):
+
+    if(color == "black"):
+
+        print("\033[30m{}\033[0m".format(text))
+
+    elif(color == "red"):
+
+        print("\033[31m{}\033[0m".format(text))
+
+    elif(color == "green"):
+
+        print("\033[32m{}\033[0m".format(text))
+
+    elif(color == "yellow"):
+
+        print("\033[33m{}\033[0m".format(text))
+
+    elif(color == "blue"):
+
+        print("\033[34m{}\033[0m".format(text))
+
+    elif(color == "magenta"):
+
+        print("\033[35m{}\033[0m".format(text))
+
+    elif(color == "cyan"):
+
+        print("\033[36m{}\033[0m".format(text))
+
+    elif(color == "white"):
+
+        print("\033[37m{}\033[0m".format(text))
+
+    else:
+
+        print("{}".format(text))
+
+
+# ---
+
+def stdout_col(text, color):
+
+    if(color == "black"):
+
+        sys.stdout.write("\033[30m{}\033[0m".format(text))
+
+    elif(color == "red"):
+
+        sys.stdout.write("\033[31m{}\033[0m".format(text))
+
+    elif(color == "green"):
+
+        sys.stdout.write("\033[32m{}\033[0m".format(text))
+
+    elif(color == "yellow"):
+
+        sys.stdout.write("\033[33m{}\033[0m".format(text))
+
+    elif(color == "blue"):
+
+        sys.stdout.write("\033[34m{}\033[0m".format(text))
+
+    elif(color == "magenta"):
+
+        sys.stdout.write("\033[35m{}\033[0m".format(text))
+
+    elif(color == "cyan"):
+
+        sys.stdout.write("\033[36m{}\033[0m".format(text))
+
+    elif(color == "white"):
+
+        sys.stdout.write("\033[37m{}\033[0m".format(text))
+
+    else:
+
+        sys.stdout.write("{}".format(text))
+
+# ---
+
+