From 7b280c004a13912fa04b7c2ba0ca59c7e9f2001e Mon Sep 17 00:00:00 2001
From: Sergey Yurchenko <sergei.yurchenko@gmail.com>
Date: Sat, 9 Sep 2023 23:32:02 +0100
Subject: [PATCH] fix the control problem where basis can overright J in
 control step 3

---
 fields.f90       | 64 ++++++++++++++++++++++++++++++++++++------------
 perturbation.f90 |  6 ++---
 2 files changed, 52 insertions(+), 18 deletions(-)

diff --git a/fields.f90 b/fields.f90
index 670d9d4..fd66a1b 100644
--- a/fields.f90
+++ b/fields.f90
@@ -1463,18 +1463,6 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
                 !
              enddo
              !
-             !call readi(i_t)
-             !
-             !do while (i_t/=0.and.i<=size(job%select_gamma))
-             !   !
-             !   i = i+1
-             !   !
-             !   job%select_gamma(i_t) = .true.
-             !   !
-             !   call readi(i_t)
-             !   !
-             !enddo 
-             !
            case ("GRAM-SCHMIDT","SVD","SCHMIDT")
              !
              job%orthogonalizer =trim(w)
@@ -2005,7 +1993,7 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
                 call readi(job%bset(imode)%range(1))
                 call readi(job%bset(imode)%range(2))
                 !
-                ! in case the range(2) is given for imode=0 we use Jrot as range(2) in order 
+                ! in case the range(2) is given for imode=0 we use Jrot as range(2)
                 !
                 if (imode==0.and.Jrot/=0) then 
                   !
@@ -2017,8 +2005,6 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
                 if (job%bset(imode)%dvrpoints==0) job%bset(imode)%dvrpoints = job%bset(imode)%range(2)+1
                 !
               case("JROT")
-                !
-                ! we use range(1) to store the Jrot value 
                 !
                 call readi(i_t)
                 !
@@ -2026,6 +2012,8 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
                   Jrot = i_t
                 endif
                 !
+                ! we use range(1) to store the Jrot value 
+                !
                 job%bset(imode)%range(1) = Jrot
                 !
               case("KROT","K")
@@ -3220,6 +3208,7 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
              case("3")
                !
                call readi(Jrot)
+               job%bset(0)%range(1) = Jrot
                !
              case("4")
                !
@@ -3228,6 +3217,51 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
                !
              end select
              !
+           case ("GAMMA")
+             !
+             if (.not.symmetry_defined) then 
+                !
+                write (out,"('FLinput: keyword GAMMA in CONTRACT cannot appear before symmetry is defined')") 
+                stop 'FLinput - symmetry should be defined before GAMMA '
+                !
+             endif
+             !
+             if (Nitems>sym%Nrepresen+1.or.Nitems==1) then  
+               write (out,"('FLinput: illegal number of irreps in gamma in DIAGONALIZER: ',i7)") Nitems-1
+               stop 'FLinput - illegal number of gammas in DIAGONALIZER'
+             endif
+             !
+             i = 0
+             job%select_gamma = .false.
+             !
+             do while (item<Nitems.and.i<size(job%select_gamma))
+                !
+                i = i + 1
+                !
+                call readu(w)
+                !
+                if (trim(w)/="-") then
+                  !
+                  read(w,*) i_t
+                  job%select_gamma(i_t) = .true.
+                  !
+                else
+                  !
+                  call readi(i_tt)
+                  !
+                  do while (i_t<i_tt.and.i<size(job%select_gamma))
+                    !
+                    i_t = i_t + 1
+                    job%select_gamma(i_t) = .true.
+                    i = i + 1
+                    !
+                  enddo
+                  i = i - 1
+                  !
+                endif
+                !
+             enddo
+             !
            case default
              !
              call report ("Unrecognized unit name "//trim(w),.true.)
diff --git a/perturbation.f90 b/perturbation.f90
index 84d4cb4..61273b8 100644
--- a/perturbation.f90
+++ b/perturbation.f90
@@ -477,7 +477,7 @@ subroutine PTinit(NPTorder,Nmodes,Npolyads)
     integer(ik)                   :: alloc,i
     integer(ik)                   :: imode,ispecies,Nspecies,Nclasses,iclass
 
-    if (job%verbose>=5) write(out,"(/'PTinit/start: Initialization of the perturbation theory elements ')") 
+    if (job%verbose>=5) write(out,"(/'PTinit/start: Initialization of the variational theory elements ')") 
 
 
     ! Define the PT-elements 
@@ -8361,10 +8361,10 @@ subroutine PThamiltonian_contract(jrot)
         !
         ! Diagonalization
         !
-        deallocate (smat(isym)%coeffs)
+        deallocate (smat(isym)%coeffs)  
         call Arraystop('PThamiltonian_contract:smat'//sym%label(isym))
         !
-        if (job%verbose>=3) write(out,"('Ready fo diagonalization...')")
+        if (job%verbose>=3) write(out,"('Ready for diagonalization...')")
         !
         if (job%verbose>=1) then 
            write (out,"(//'Size of the symmetrized hamiltonian = ',i7,' Symmetry = ',a4)") dimen_s,sym%label(isym)