WIP: Attempt for better initialize GLE

Does not seem to work so commented out for now.
Trying to solve #25

abin.F90

@@ -100,7 +100,7 @@ program abin_dyn
 
 !  End of transformations
 
-!-----Note that amt equals am if staging is off
+   ! Note that 'amt' equals 'am' for non-PI simulations
    px = amt * vx
    py = amt * vy
    pz = amt * vz

force_tera.F90

@@ -404,14 +404,14 @@ subroutine connect_terachem( itera )
 
    end if
 
-   write(6,'(2a)') 'Establishing connection to TeraChem port: ', trim(port_name)
+   write(6,'(a)') 'Establishing connection to TeraChem...'
    ! ----------------------------------------
    ! Establish new communicator via port name
    ! ----------------------------------------
    call flush(6)
    call MPI_COMM_CONNECT(port_name, MPI_INFO_NULL, 0, MPI_COMM_SELF, newcomm, ierr)
    call handle_mpi_error(ierr)
-   write(6,'(a,i0)') 'Established new communicator:', newcomm
+   write(6,'(a,i0)') 'Established a new communicator:', newcomm
 
    newcomms(itera) = newcomm
 

gle.F90

@@ -243,18 +243,19 @@ subroutine gle_init(dt)
    end if
    close(122)
 
+   ! WARNING: gA is rewritten here
+   ! TODO: do not overwrite gA
    call compute_propagator(gA, gC, gT, gS, dt)
 
-
-   ! then, we must initialize the auxiliary vectors. we keep general - as we might be 
-   ! using non-diagonal C to break detailed balance - and we use cholesky decomposition
-   ! of C. again, since one would like to initialize correctly the velocities in 
-   ! case of generic C, we use an extra slot for gp for the physical momentum, as we 
-   ! could then use it to initialize the momentum in the calling code
+   ! Initialize the auxiliary vectors. 
+   ! we keep general - as we might be using non-diagonal C
+   ! to break detailed balance - and we use cholesky decomposition of C
+   ! since one would like to initialize correctly the velocities in 
 
    ! DH: ps rewritten in init.f90 if irest.eq.1
-   ! always initialize, maybe rewritten in restart
 
+   ! TODO: Do not overwrite gA
+   ! TODO: Move this inside initialize_momenta
    gA = gC   
    call cholesky(gA, gC, ns+1)
 
@@ -286,9 +287,9 @@ subroutine gle_init(dt)
    end subroutine gle_init
 
    subroutine initialize_momenta(C, iw)
-   !use mod_arrays,  only: px, py, pz
+   !use mod_arrays,  only: px, py, pz, vx, vy, vz, amt
    use mod_general, only: natom
-   use mod_utils,   only: abinerror
+   use mod_utils,   only: abinerror, print_xyz_arrays
    real(DP), intent(in) :: C(:,:)
    integer, intent(in)  :: iw
    real(DP),allocatable :: gr(:)
@@ -301,14 +302,27 @@ subroutine initialize_momenta(C, iw)
 
    allocate(gr(ns+1))
 
-   do j=1,natom*3
+   do j = 1, natom * 3
       call gautrg(gr, ns + 1)
-      ! TODO, actually pass this to initialize momenta
       gp(j,:) = matmul(C, gr)
    end do
 
-   do j=1,natom*3
-      do i=1,ns
+   ! TODO, actually pass this to initialize momenta
+   ! case of generic C, we use an extra slot for gp for the physical momentum
+   !do i = 1, natom
+   !   px(i,iw) = gp(i,1)           
+   !   py(i,iw) = gp(i + natom,1)     
+   !   pz(i,iw) = gp(i + 2*natom,1)   
+   !end do
+   !vx = px / amt
+   !vy = py / amt
+   !vz = pz / amt
+
+   !call print_xyz_arrays(px, py, pz)
+   !call print_xyz_arrays(vx, vy, vz)
+
+   do j = 1, natom*3
+      do i = 1, ns
          ps(j,i,iw) = gp(j,i+1)
       enddo
    enddo
@@ -330,17 +344,18 @@ subroutine finalize_gle()
    end subroutine finalize_gle
 
 
+   ! Matrix A is rewritten on output
    subroutine compute_propagator(A, C, T, S, dt)
-      ! Matrix A is rewritten on output
    real(DP), intent(inout) :: A(:,:) 
    real(DP), intent(in)    :: C(:,:)
    real(DP), intent(out)   :: T(:,:), S(:,:)
    real(DP), intent(in)    :: dt
 
-   ! the deterministic part of the propagator is obtained in a second
+   ! the deterministic part of the propagator
    call matrix_exp(-dt*A, ns+1, 15, 15, T)
 
-   ! the stochastic part is just as easy. we use gA as a temporary array
+   ! the stochastic part, we use A as a temporary array
+   ! TODO: Do not overwrite A, makes things confusing
    A = C - matmul(T, matmul(C, transpose(T)) )
    call cholesky(A, S, ns+1)
 
@@ -477,8 +492,8 @@ subroutine matrix_exp(M, n, j, k, EM)
     enddo
   end subroutine matrix_exp
 
-  ! TODO: replace by more stable procedure from i-Py???
-  ! brute-force "stabilized" cholesky decomposition.
+  ! TODO: replace by more stable procedure from i-Py.
+  ! Brute-force "stabilized" cholesky decomposition.
   ! in practice, we compute LDL^T decomposition, and force
   ! to zero negative eigenvalues.
   subroutine cholesky(SST, S, n)

init.F90

@@ -62,7 +62,7 @@ subroutine init(dt)
    character(len=60)   :: chdivider
    character(len=60)   :: mdtype
    LOGICAL :: file_exists
-   logical        :: rem_comvel, rem_comrot
+   logical :: rem_comvel, rem_comrot
 !  real(DP) :: wnw=5.0d-5
    ! Used for MPI calls
    integer :: ierr
@@ -505,11 +505,11 @@ subroutine init(dt)
 
       if (my_rank.eq.0)then
          if (temp0.gt.0)then
-            write(*,*)'Initial temperature in Kelvins =', temp0
+            write(*,*)'Initial temperature [K] =', temp0
          else
-            write(*,*)'Initial temperature in Kelvins =', temp
+            write(*,*)'Initial temperature [K] =', temp
          end if
-         if (inose.ne.0) write(*,*)'Target temperature in Kelvins =', temp
+         if (inose.ne.0) write(*,*)'Target temperature [K] =', temp
       end if
 
       ! conversion of temperature from K to au
@@ -574,10 +574,11 @@ subroutine init(dt)
 
 !     SETTING initial velocities according to the Maxwell-Boltzmann distribution
       if(irest.eq.0.and.chveloc.eq.'')then
+         ! TODO: GLE thermostat, initialize momenta in gle_init
          if (temp0.ge.0)then
-            call vinit(TEMP0, am, vx, vy, vz)
+            call vinit(temp0, am, vx, vy, vz)
          else
-            call vinit(TEMP, am, vx, vy, vz)
+            call vinit(temp, am, vx, vy, vz)
          end if
       end if
 
@@ -624,12 +625,11 @@ subroutine init(dt)
       ! Otherwise, just print the temperature.
       call ScaleVelocities(vx, vy, vz)
 
-
-!-----some stuff for spherical boundary onditions
+      ! Initialize spherical boundary onditions
       if(isbc.eq.1) call sbc_init(x,y,z)
 
-!-----inames initialization for the MM part. 
-!-----We do this also because string comparison is very costly
+      ! inames initialization for the MM part. 
+      ! We do this also because string comparison is very costly
       if(iqmmm.eq.3.or.pot.eq.'mm') allocate( inames(natom) )
 
       if(iqmmm.eq.3.or.pot.eq.'mm')then