Merge pull request #6 from lgcrego/alpha-SO

adapting dynamics routines to spin handling and output
lgcrego · Dec 19, 2020 · 17b1e21 · 17b1e21
2 parents 1d457f6 + cf5b05b
commit 17b1e21
Show file tree

Hide file tree

Showing 12 changed files with 321 additions and 275 deletions.
diff --git a/AO_adiabatic.f b/AO_adiabatic.f
@@ -16,7 +16,7 @@ module AO_adiabatic_m
                                              DensityMatrix, AutoCorrelation,  &
                                              CT_dump_step, Environ_step,      &
                                              driver, HFP_Forces ,             &
-                                             step_security
+                                             step_security , SOC
     use Babel_m                     , only : Coords_from_Universe, trj, MD_dt                            
     use Allocation_m                , only : Allocate_UnitCell ,              &
                                              DeAllocate_UnitCell ,            &
@@ -38,7 +38,8 @@ module AO_adiabatic_m
     use Schroedinger_m              , only : DeAllocate_QDyn
     use Psi_Squared_Cube_Format     , only : Gaussian_Cube_Format
     use Data_Output                 , only : Populations ,                    &
-                                             Net_Charge                       
+                                             Net_Charge ,                     &
+                                             FileName                      
     use Backup_m                    , only : Security_Copy ,                  &
                                              Restart_state ,                  &
                                              Restart_Sys                      
@@ -123,7 +124,7 @@ subroutine AO_adiabatic( Qdyn , final_it )
     CALL DUAL_wvpckts
 
     ! save populations(t + t_rate)  and  update Net_Charge ...
-    QDyn%dyn(it,:,:) = Populations( QDyn%fragments , ExCell_basis , DUAL_bra , DUAL_ket , t )
+    QDyn%dyn(it,:,:,:) = Populations( QDyn%fragments , ExCell_basis , DUAL_bra , DUAL_ket , t )
 
     if( mod(it,CT_dump_step) == 0 ) CALL dump_Qdyn( Qdyn )
 
@@ -303,7 +304,7 @@ subroutine Preprocess( QDyn )
 CALL DUAL_wvpckts
 
 ! save populations ...
-QDyn%dyn(it,:,:) = Populations( QDyn%fragments , ExCell_basis , DUAL_bra , DUAL_ket , t_i )
+QDyn%dyn(it,:,:,:) = Populations( QDyn%fragments , ExCell_basis , DUAL_bra , DUAL_ket , t_i )
 CALL dump_Qdyn( Qdyn )
 
 If( GaussianCube ) CALL Send_to_GaussianCube( it )
@@ -501,40 +502,45 @@ subroutine dump_Qdyn( Qdyn )
 type(f_time)    , intent(in) :: QDyn
 
 ! local variables ...
-integer    :: nf , n
-complex*16 :: wp_energy(n_part)
-
-do n = 1 , n_part
-
-    if( eh_tag(n) == "XX" ) cycle
-
-    wp_energy(n) = sum(MO_bra(:,n)*UNI%erg(:)*MO_ket(:,n)) 
-
-    If( it == 1 ) then
-
-        open( unit = 52 , file = "dyn.trunk/"//eh_tag(n)//"_survival.dat" , status = "replace" , action = "write" , position = "append" )
-        write(52,11) "#" ,( nf+1 , nf=0,size(QDyn%fragments)+1 )  ! <== numbered columns for your eyes only ...
-        write(52,12) "#" , QDyn%fragments , "total"
-
-        open( unit = 53 , file = "dyn.trunk/"//eh_tag(n)//"_wp_energy.dat" , status = "replace" , action = "write" , position = "append" )
-
-    else
-
-        open( unit = 52 , file = "dyn.trunk/"//eh_tag(n)//"_survival.dat"  , status = "unknown", action = "write" , position = "append" )
-        open( unit = 53 , file = "dyn.trunk/"//eh_tag(n)//"_wp_energy.dat" , status = "unknown", action = "write" , position = "append" )
-
-    end If
-
-    ! dumps el-&-hl populations ...
-    write(52,13) ( QDyn%dyn(it,nf,n) , nf=0,size(QDyn%fragments)+1 ) 
-
-    ! dumps el-&-hl wavepachet energies ...
-    write(53,14) QDyn%dyn(it,0,n) , real( wp_energy(n) ) , dimag( wp_energy(n) )
-
-    close(52)
-    close(53)
-
-end do
+integer                       :: nf , n , spin , n_spin
+complex*16                    :: wp_energy(n_part)
+character(len=:) ,allocatable :: f_name
+
+n_spin = merge(2,1,SOC)
+
+do spin = 1 , n_spin
+   do n = 1 , n_part
+
+         if( eh_tag(n) == "XX" ) cycle
+
+         wp_energy(n) = sum(MO_bra(:,n)*UNI%erg(:)*MO_ket(:,n)) 
+
+         If( it == 1 ) then
+             call FileName( f_name , n , spin , instance="dens" )
+             open( unit = 52 , file = f_name , status = "replace" , action = "write" , position = "append" )
+             write(52,11) "#" ,( nf+1 , nf=0,size(QDyn%fragments)+1 )  ! <== numbered columns for your eyes only ...
+             write(52,12) "#" , QDyn%fragments , "total"
+
+             call FileName( f_name , n , spin , instance="erg" )
+             open( unit = 53 , file = f_name , status = "replace" , action = "write" , position = "append" )
+         else   
+             call FileName( f_name , n , spin , instance="dens" )
+             open( unit = 52 , file = f_name , status = "unknown" , action = "write" , position = "append" )
+             call FileName( f_name , n , spin , instance="erg" )
+             open( unit = 53 , file = f_name , status = "unknown" , action = "write" , position = "append" )
+         end If
+
+         ! dumps el-&-hl populations ...
+         write(52,13) ( QDyn%dyn(it,nf,n,spin) , nf=0,size(QDyn%fragments)+1 )
+
+         ! dumps el-&-hl wavepachet energies ...
+         write(53,14) QDyn%dyn(it,0,n,spin) , real( wp_energy(n) ) , dimag( wp_energy(n) )         
+
+         close(52)
+         close(53)
+
+         end do
+         end do
 
 ! QM_erg = E_occ - E_empty ; to be used in MM_dynamics energy balance ...
 Unit_Cell% QM_erg = real( wp_energy(1) ) - real( wp_energy(2) )

diff --git a/ElHl_Chebyshev.f b/ElHl_Chebyshev.f
@@ -137,7 +137,7 @@ subroutine preprocess_ElHl_Chebyshev( system , basis , AO_bra , AO_ket , Dual_br
 !==============================================
 
 ! save populations(time=t_i) ...
-QDyn%dyn(it,:,:) = Populations( QDyn%fragments , basis , DUAL_bra , DUAL_ket , t_i )
+QDyn%dyn(it,:,:,:) = Populations( QDyn%fragments , basis , DUAL_bra , DUAL_ket , t_i )
 CALL dump_Qdyn( Qdyn , it )
 
 !==============================================
@@ -272,7 +272,7 @@ subroutine ElHl_Chebyshev( system , basis , AO_bra , AO_ket , Dual_bra , Dual_ke
 CALL store_Hprime( N , H_prime ) 
 
 ! save populations(time) ...
-QDyn%dyn(it,:,:) = Populations( QDyn%fragments , basis , DUAL_bra , DUAL_ket , t )
+QDyn%dyn(it,:,:,:) = Populations( QDyn%fragments , basis , DUAL_bra , DUAL_ket , t )
 
 if( mod(it,CT_dump_step) == 0 ) CALL dump_Qdyn( Qdyn , it )
 

diff --git a/ElHl_schroedinger.f b/ElHl_schroedinger.f
@@ -7,13 +7,12 @@ module Schroedinger_m
  use parameters_m               , only : t_i , t_f , n_t , n_part , GaussianCube , preview, &
                                          GaussianCube_step ,  DP_Moment , electron_state ,  &
                                          Coulomb_ , DensityMatrix , driver , SOC
- use tuning_m                   , only : spin_tag 
  use Allocation_m               , only : Allocate_Brackets , DeAllocate_Structures
  use Babel_m                    , only : trj , Coords_from_Universe
  use Structure_Builder          , only : Unit_Cell , Extended_Cell , Generate_Structure
  use FMO_m                      , only : FMO_analysis , orbital , eh_tag
  use DP_main_m                  , only : Dipole_Moment
- use Data_Output                , only : Populations , Net_Charge
+ use Data_Output                , only : Populations , Net_Charge , FileName
  use Psi_Squared_Cube_Format    , only : Gaussian_Cube_Format
  use Backup_m                   , only : Security_Copy , Restart_state
  use Auto_Correlation_m         , only : MO_Occupation
@@ -24,10 +23,11 @@ module Schroedinger_m
     private
 
     ! module variables ...
-    Complex*16   , ALLOCATABLE , dimension(:,:)   :: MO_bra , MO_ket , AO_bra , AO_ket , DUAL_ket , DUAL_bra
-    Real*8       , ALLOCATABLE , dimension(:,:,:) :: Pops
-    type(f_time)                                  :: Mean_QDyn , aux
-    integer                                       :: iter = 0 
+    Complex*16   , ALLOCATABLE , dimension(:,:)     :: MO_bra , MO_ket , AO_bra , AO_ket , DUAL_ket , DUAL_bra
+    Real*8       , ALLOCATABLE , dimension(:,:,:,:) :: Pops
+    type(f_time)                                    :: Mean_QDyn , aux
+    integer                                         :: n_spin
+    integer                                         :: iter = 0 
 
  contains
 !
@@ -51,7 +51,7 @@ subroutine Simple_dynamics(system, basis, UNI, QDyn )
 
 ! ------------------ preprocess stuff --------------------
 
-allocate( Pops( n_t , 0:size(system%list_of_fragments)+1 , n_part ) )
+allocate( Pops(n_t , 0:size(system%list_of_fragments)+1 , n_part , n_spin) )
 
 mm = size(basis) ; nn = n_part
 
@@ -98,9 +98,9 @@ subroutine Simple_dynamics(system, basis, UNI, QDyn )
 CALL DZgemm( 'T' , 'N' , mm , nn , mm , C_one , UNI%L , mm , MO_bra , mm , C_zero , DUAL_bra , mm )
 
 ! save populations ...
-Pops(1,:,:) = Populations( QDyn%fragments , basis , DUAL_bra , DUAL_ket , t_i )
+Pops(1,:,:,:) = Populations( QDyn%fragments , basis , DUAL_bra , DUAL_ket , t_i )
 
-QDyn%dyn(1,:,:) = Pops(1,:,:)
+QDyn%dyn(1,:,:,:) = Pops(1,:,:,:)
 
 If( DensityMatrix ) then
     If( n_part == 1 ) CALL MO_Occupation( t_i, MO_bra, MO_ket, UNI )
@@ -143,8 +143,8 @@ subroutine Simple_dynamics(system, basis, UNI, QDyn )
     CALL dzgemm( 'T' , 'N' , mm , nn , mm , C_one , UNI%L , mm , MO_bra , mm , C_zero , DUAL_bra , mm )
 
     ! get populations ...
-    Pops(it,:,:) = Populations( QDyn%fragments , basis , DUAL_bra , DUAL_ket , t )
-    QDyn%dyn(it,:,:) = Pops(it,:,:)
+    Pops(it,:,:,:) = Populations( QDyn%fragments , basis , DUAL_bra , DUAL_ket , t )
+    QDyn%dyn(it,:,:,:) = Pops(it,:,:,:)
 
     ! LOCAL representation for film STO production ...
     AO_bra = DUAL_bra
@@ -213,7 +213,7 @@ subroutine RunningStat( Qdyn , instance )
 If( present(instance) ) then
    Qdyn% std = sqrt( Qdyn% std/float(iter-1) )
    ! preserving time; does not undergo std ...
-   Qdyn% std(:,0,:) = Qdyn% dyn(:,0,:)  
+   Qdyn% std(:,0,:,:) = Qdyn% dyn(:,0,:,:)  
 
    return
 end If
@@ -253,12 +253,10 @@ subroutine dump_Qdyn( Qdyn , UNI )
 type(R_eigen)   , intent(in) :: UNI
 
 ! local variables ...
-integer                       :: nf , n , it , n_spin , spin
+integer                       :: nf , n , it , spin
 complex*16                    :: wp_energy
 character(len=:) ,allocatable :: f_name
 
-n_spin = merge(2,1,SOC)
-
 do spin = 1 , n_spin
 
    do n = 1 , n_part
@@ -277,10 +275,10 @@ subroutine dump_Qdyn( Qdyn , UNI )
          DO it = 1 , n_t
 
              ! dumps el-&-hl populations ...
-             write(52,13) ( QDyn%dyn(it,nf,n) , nf=0,size(QDyn%fragments)+1 )
+             write(52,13) ( QDyn%dyn(it,nf,n,spin) , nf=0,size(QDyn%fragments)+1 )
 
              ! dumps el-&-hl wavepachet energies ...
-             write(53,14) QDyn%dyn(it,0,n) , real(wp_energy) , dimag(wp_energy)
+             write(53,14) QDyn%dyn(it,0,n,spin) , real(wp_energy) , dimag(wp_energy)
 
          end do
 
@@ -308,10 +306,12 @@ subroutine DeAllocate_QDyn( QDyn , flag )
 character(*)  , intent(in)    :: flag
 
 ! local variable ...
-integer      :: N_of_fragments
+integer      :: N_of_fragments 
 character(1) :: first_in_line
 logical      :: A_flag
 
+n_spin = merge(2,1,SOC)
+
 select case( flag )
 
     case( "alloc" )
@@ -340,10 +340,10 @@ subroutine DeAllocate_QDyn( QDyn , flag )
         Extended_Cell%list_of_fragments(1) = merge( "A" , "D" , A_flag )
 
         ! QDyn%dyn = ( time ; fragments ; all fragments ) ...
-        allocate( QDyn%fragments( size(Extended_Cell % list_of_fragments) ) , source = Extended_Cell % list_of_fragments )
-        allocate( QDyn%dyn      ( n_t , 0:N_of_fragments+1 , n_part       ) , source = 0.d0                              )
+        allocate( QDyn%fragments( size(Extended_Cell % list_of_fragments)    ) , source = Extended_Cell % list_of_fragments )
+        allocate( QDyn%dyn      ( n_t , 0:N_of_fragments+1 , n_part , n_spin ) , source = 0.d0                              )
         If( driver == "avrg_confgs" ) then
-              allocate( QDyn%std( n_t , 0:N_of_fragments+1 , n_part       ) , source = 0.d0                              )
+              allocate( QDyn%std( n_t , 0:N_of_fragments+1 , n_part , n_spin ) , source = 0.d0                              )
         End If
 
         ! allocatating Net_Charte for future use ...
@@ -364,48 +364,4 @@ end subroutine DeAllocate_QDyn
 !
 !
 !
-!===============================================
-subroutine FileName( f_name , n , s , instance )
-!===============================================
-implicit none
-character(len=:) , allocatable , intent(out) :: f_name
-integer                        , intent(in)  :: n
-integer                        , intent(in)  :: s
-character(len=*)               , intent(in)  :: instance
-
-! s stands for spin flag ...
-! n stands for n_part flag ...
-
-select case ( instance )
-
-       case ("dens")
-            if( spin_tag(s) == "XX" ) then
-                allocate( character(len=21) :: f_name )
-                f_name = "dyn.trunk/"//eh_tag(n)//"_dens.dat" 
-            elseif( spin_tag(s) == "up" ) then
-                allocate( character(len=24) :: f_name )
-                f_name = "dyn.trunk/"//eh_tag(n)//"_dens_up.dat" 
-            elseif( spin_tag(s) == "dw" ) then
-                allocate( character(len=26) :: f_name )
-                f_name = "dyn.trunk/"//eh_tag(n)//"_dens_down.dat" 
-            end If
-
-       case ("erg") 
-            if( spin_tag(s) == "XX" ) then
-                allocate( character(len=23) :: f_name )
-                f_name = "dyn.trunk/"//eh_tag(n)//"_wp-erg.dat" 
-            elseif( spin_tag(s) == "up" ) then
-                allocate( character(len=26) :: f_name )
-                f_name = "dyn.trunk/"//eh_tag(n)//"_wp-erg_up.dat" 
-            elseif( spin_tag(s) == "dw" ) then
-                allocate( character(len=28) :: f_name )
-                f_name = "dyn.trunk/"//eh_tag(n)//"_wp-erg_down.dat" 
-            end If
-
-end select
-
-end subroutine FileName
-!
-!
-!
 end module Schroedinger_m
diff --git a/Taylor.f b/Taylor.f
@@ -5,11 +5,11 @@ module Taylor_m
     use blas95
     use lapack95
     use ifport
-    use parameters_m        , only : t_i, frame_step, n_part, restart, CT_dump_step                  
+    use parameters_m        , only : t_i, frame_step, n_part, restart, CT_dump_step , SOC
     use Structure_Builder   , only : Unit_Cell                                      
     use Overlap_Builder     , only : Overlap_Matrix
     use FMO_m               , only : eh_tag                  
-    use Data_Output         , only : Populations
+    use Data_Output         , only : Populations , FileName
     use Matrix_Math
 
     public  :: Propagation, dump_Qdyn
@@ -248,34 +248,42 @@ subroutine dump_Qdyn( Qdyn , it )
 integer         , intent(in) :: it 
 
 ! local variables ...
-integer :: nf , n
-
-do n = 1 , n_part
-
-    If( eh_tag(n) == "XX" ) cycle
-
-    If( it == 1 ) then
-        open( unit = 52 , file = "dyn.trunk/"//eh_tag(n)//"_survival.dat" , status = "replace" , action = "write" , position = "append" )
-        write(52,15) "#" ,( nf+1 , nf=0,size(QDyn%fragments)+1 )  ! <== numbered columns for your eyes only ...
-        write(52,12) "#" , QDyn%fragments , "total"
-
-        open( unit = 53 , file = "dyn.trunk/"//eh_tag(n)//"_wp_energy.dat" , status = "replace" , action = "write" , position = "append" )
-
-    else
-        open( unit = 52 , file = "dyn.trunk/"//eh_tag(n)//"_survival.dat"  , status = "unknown", action = "write" , position = "append" )
-        open( unit = 53 , file = "dyn.trunk/"//eh_tag(n)//"_wp_energy.dat" , status = "unknown", action = "write" , position = "append" )        
-    end If
-
-    ! dumps el-&-hl populations ...    
-    write(52,13) ( QDyn%dyn(it,nf,n) , nf=0,size(QDyn%fragments)+1 )
-
-    ! dumps el-&-hl wavepachet energies ...
-    write(53,14) QDyn%dyn(it,0,n) , real( Unit_Cell% QM_wp_erg(n) ) , dimag( Unit_Cell% QM_wp_erg(n) )
-
-    close(52)
-    close(53)
-
-end do
+integer                       :: nf , n , spin , n_spin
+character(len=:) ,allocatable :: f_name
+
+n_spin = merge(2,1,SOC)
+
+do spin = 1 , n_spin
+   do n = 1 , n_part
+
+         if( eh_tag(n) == "XX" ) cycle
+
+         If( it == 1 ) then
+             call FileName( f_name , n , spin , instance="dens" )
+             open( unit = 52 , file = f_name , status = "replace" , action = "write" , position = "append" )
+             write(52,15) "#" ,( nf+1 , nf=0,size(QDyn%fragments)+1 )  ! <== numbered columns for your eyes only ...
+             write(52,12) "#" , QDyn%fragments , "total"
+
+             call FileName( f_name , n , spin , instance="erg" )
+             open( unit = 53 , file = f_name , status = "replace" , action = "write" , position = "append" )
+         else   
+             call FileName( f_name , n , spin , instance="dens" )
+             open( unit = 52 , file = f_name , status = "unknown" , action = "write" , position = "append" )
+             call FileName( f_name , n , spin , instance="erg" )
+             open( unit = 53 , file = f_name , status = "unknown" , action = "write" , position = "append" )
+         end If
+
+         ! dumps el-&-hl populations ...
+         write(52,13) ( QDyn%dyn(it,nf,n,spin) , nf=0,size(QDyn%fragments)+1 )
+
+         ! dumps el-&-hl wavepachet energies ...
+         write(53,14) QDyn%dyn(it,0,n,spin) , real( Unit_Cell% QM_wp_erg(n) ) , dimag( Unit_Cell% QM_wp_erg(n) )
+
+         close(52)
+         close(53)
+
+         end do
+         end do
 
 12 FORMAT(/15A10)
 13 FORMAT(F11.6,14F10.5)