CSNV, bug in contracted basis for sing_triatom, compact form for fiel…

…ds and energies, fix, fit keywords, H2CO PES with isotopologues

fields.f90

@@ -665,7 +665,7 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
    logical :: eof,zmat_defined,basis_defined,equil_defined,pot_defined,symmetry_defined,extF_defined,refer_defined,chk_defined
    logical :: kinetic_defined,pot_form_compact = .false.,extF_form_compact = .false.
    character(len=cl) :: Molecule,pot_coeff_type,exfF_coeff_type,chk_type,controlstep
-   character(len=wl) :: w,ioname
+   character(len=wl) :: w,ioname,w_t
    real(rk)    :: lfact,f_t, func_coef
    integer(ik) :: i,iatom,imode, ifunc,numterms,  numfunc, in_expo, out_expo ,natoms,alloc,Nparam,iparam,i_t,i_tt
    integer(ik) :: Nbonds,Nangles,Ndihedrals,j,ispecies,imu,iterm,Ncoords,icoords
@@ -681,7 +681,7 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
    logical           :: zEchoInput = .true.
    integer,parameter :: max_input_lines=500000  ! maximum length (in lines) of input. 500,000 lines is plenty..
    !                                              ! to avoid feeding in GB of data by mistake.
-   integer(ik)  :: Nparam_check    !number of parameters as determined automatically by duo (Nparam is specified in input).
+   !integer(ik)  :: Nparam_check    !number of parameters as determined automatically by duo (Nparam is specified in input).
    !
    !
    ! default values: 
@@ -3565,40 +3565,7 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
          !
          ! find the number of points/parameters in input
          !
-         Nparam_check = 0
-         !
-         call input_options(echo_lines=.false.)
-         !
-         do while (trim(w)/="END")
-            !
-            call read_line(eof,iut) ; if (eof) exit
-            !
-            call readu(w)
-            !
-            Nparam_check = Nparam_check+1
-            !
-         enddo
-         !
-         Nparam_check = Nparam_check-1
-         !
-         if (zEchoInput) call input_options(echo_lines=.true.)
-         !
-         if (trim(w) /= "END") then
-             call report ("ERROR: Cannot find `END' statement)",.true.)
-         endif
-         !
-         ! go back to beginning of VALUES block and reset `w' to original value
-         do i=1, Nparam_check+1
-           backspace(unit=iut)
-         enddo
-         !
-         w = trim(pot_coeff_type)
-         !
-         Nparam = Nparam_check
-         !
-         if (Nparam <= 0) then
-             call report ("ERROR: Number of points or parameters <= 0 )",.true.)
-         endif
+         call count_parameters_in_input_file(iut,Nparam)
          !
          ! Allocation of the pot. parameters 
          !
@@ -3714,6 +3681,8 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
             !
          enddo
          !
+         call count_parameters_in_input_file(iut,Nparam)
+         !
          ! Allocation of the pot. parameters 
          !
          allocate (molec%mep_params(Nparam),stat=alloc)
@@ -4289,9 +4258,9 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
              !
              call readf(fitting%threshold_coeff)
              !
-           case('OBS','OBS_ENERGIES')
+           case('OBS','OBS_ENERGIES','ENERGIES')
              !
-             call readi(fitting%Nenergies)
+             call count_parameters_in_input_file(iut,fitting%Nenergies)
              !
              allocate (fitting%obs(1:fitting%Nenergies),stat=alloc)
              if (alloc/=0) then
@@ -4625,6 +4594,13 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
             !
          enddo
          !
+         call count_parameters_in_input_file(iut,Nparam)
+         !
+         if ( Nparam /= sum(extF%nterms(:)) ) then
+             write (out,"('FLinput: number of rows in EXTERNAL <> NPARAM',i8,' <> sum of ',20i5)") Nparam,extF%nterms(:)
+             stop 'FLinput - illigal number of rows (parameters) in ExtF'
+         endif 
+         !
          ! Allocation of the pot. parameters 
          !
          allocate (extF%coef(Nparam,extF%rank),extF%name(Nparam,extF%rank),& 
@@ -4639,9 +4615,7 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
          extF%coef = 0
          extF%name = 'xxxx'
          extF%term = -1
-         extF%ifit = 1
-         !
-         Nparam = sum(extF%nterms(:))
+         extF%ifit = 0
          !
          do imu = 1,extF%rank
            !
@@ -4667,6 +4641,29 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
                !
                call readf(f_t) ; extF%coef(iterm,imu) = f_t
                !
+               if (extF_form_compact.and.nitems>2) then
+                  !
+                  call readu(w_t)
+                  !
+                  select case (trim(w_t))
+                    !
+                  case("FIT")
+                    !
+                    ! positv index means the parameter is fitted 
+                    !
+                    extF%ifit(iterm,imu) = 1
+                    !
+                  case("FIX")
+                    !
+                    ! negative index means the value of the paramater is fixed to the input value
+                    ! zero is for non fitted parameters and non-fixed values
+                    !
+                    extF%ifit(iterm,imu) = -1
+                    !
+                  end select
+                  !
+               endif
+               !
              case("POWERS")
                !
                if ((extF_form_compact.and.nitems<trove%Ncoords+2).or.&
@@ -4686,6 +4683,29 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
                if (.not.extF_form_compact) call readi(extF%ifit(iterm,imu))
                call readf(f_t); extF%coef(iterm,imu) = f_t
                !
+               if (extF_form_compact.and.nitems>trove%Ncoords+2) then
+                  !
+                  call readu(w_t)
+                  !
+                  select case (trim(w_t))
+                    !
+                  case("FIT")
+                    !
+                    ! positv index means the parameter is fitted 
+                    !
+                    extF%ifit(iterm,imu) = 1
+                    !
+                  case("FIX")
+                    !
+                    ! negative index means the value of the paramater is fixed to the input value
+                    ! zero is for non fitted parameters and non-fixed values
+                    !
+                    extF%ifit(iterm,imu) = -1
+                    !
+                  end select
+                  !
+               endif
+               !
                write(my_fmt,'(a,i0,a)') "(a,",Ncoords,"i1)"
                !
                write(extF%name(iterm,imu),my_fmt) 'f',(mod(extF%term(i,iterm,imu),10),i=1,trove%Ncoords)
@@ -4932,6 +4952,50 @@ subroutine print_symmetries
      !
    end subroutine print_symmetries
    !
+   subroutine count_parameters_in_input_file(iut,Nparam_check)
+     !
+     integer(ik),intent(in)  :: iut
+     integer(ik),intent(out) :: Nparam_check
+     logical :: eof
+     character(len=wl) :: w
+     integer(ik) :: i 
+         !
+         Nparam_check = 0
+         !
+         w = ""
+         !
+         call input_options(echo_lines=.false.)
+         !
+         do while (trim(w)/="END")
+            !
+            call read_line(eof,iut) ; if (eof) exit
+            !
+            call readu(w)
+            !
+            Nparam_check = Nparam_check+1
+            !
+         enddo
+         !
+         Nparam_check = Nparam_check-1
+         !
+         if (Nparam_check <= 0) then
+             call report ("ERROR: Number of parameters or energies <= 0 )",.true.)
+         endif
+         !
+         if (zEchoInput) call input_options(echo_lines=.true.)
+         !
+         if (trim(w) /= "END") then
+             call report ("ERROR: Cannot find `END' statement)",.true.)
+         endif
+         !
+         ! go back to beginning of VALUES block and reset `w' to original value
+         do i=1, Nparam_check+1
+           backspace(unit=iut)
+         enddo
+         !
+      end subroutine count_parameters_in_input_file
+
+   !
 end subroutine FLReadInput
 
 subroutine check_read_save_none(w,place)
@@ -18632,7 +18696,7 @@ subroutine FLbset1DNew(ibs,BSsize)
                 stop 'FLbset1DNew, phivphi_t  - out of memory'
              end if
              !
-             !$omp do private(vl,nl,krot1,k_l,vr,nr,krot2,k_r,Tcoeff,iterm,k1,k2,imu,mat_t) schedule(static)
+             !$omp do private(vl,k,nl,krot1,k_l,vr,nr,krot2,k_r,Tcoeff,iterm,k1,k2,imu,mat_t) schedule(static)
              do vl = 0,bs%Size
                 !
                 read (io_slot,rec=vl+1) (phil_leg(k),k=0,npoints),(dphil_leg(k),k=0,npoints)

kin_xy2.f90

@@ -433,8 +433,8 @@ subroutine MLkinetic_xyz_bisect_EKE(nmodes,Nterms,rho,g_vib,g_rot,g_cor,pseudo)
    real(ark),parameter  :: rho_threshold = 0.02_rk
      !
      if (manifold/=1) then
-       write(out,"('MLkinetic_xyz_bond_EKE-error: can be used with non-rigid case only')")
-       stop 'MLkinetic_xyz_bond_EKE can be used only with npoints>0'
+       write(out,"('MLkinetic_xyz_bisect_EKE-error: can be used with non-rigid case only')")
+       stop 'MLkinetic_xyz_bisect_EKE can be used only with npoints>0'
      endif
      !
      mX = molec%AtomMasses(1)

mol_xy2.f90

@@ -2765,7 +2765,7 @@ subroutine ML_symmetry_transformation_XY2(ioper,natoms,src,dst)
           write (out,"('ML_symmetry_transformation_XY2: symmetry ',a,' unknown')") trim(molec%symmetry)
           stop 'ML_symmetry_transformation_XY2 - bad symm. type'
           !
-       case('CS','CS(M)')
+       case('CS','CS(M)','CSN','CSN(M)')
           !
           dst = src
           !
@@ -2841,6 +2841,23 @@ subroutine ML_rotsymmetry_XY2(J,K,tau,gamma,ideg)
            !
          endif 
          !
+      case('CSN')
+         !
+         gamma = 0 
+         ideg = 1
+         !
+         if (molec%AtomMasses(2)/=molec%AtomMasses(3).or.molec%req(1)/=molec%req(2)) then
+           !
+           if (mod(tau+2,2)==0) gamma = 1+2*K !; return
+           if (mod(tau+2,2)/=0) gamma = 2+2*K !; return
+           !
+         else
+           !
+           if (mod(K+tau+2,2)==0) gamma = 1+2*K  !; return
+           if (mod(K+tau+2,2)/=0) gamma = 2+2*K  !; return
+           !
+         endif 
+         !
       case('D2H(M)')
          !
          gamma = 0 
@@ -3093,6 +3110,23 @@ subroutine ML_rotsymmetry_XY2(J,K,tau,gamma,ideg)
          if (mod(K+2,2)/=0.and.tau==0) gamma = 2+4*K
          if (mod(K+2,2)/=0.and.tau==1) gamma = 3+4*K
          !
+      case('CSN')
+         !
+         gamma = 0 
+         ideg = 1
+         !
+         if (molec%AtomMasses(2)/=molec%AtomMasses(3).or.molec%req(1)/=molec%req(2)) then
+           !
+           if (mod(tau+2,2)==0) gamma = 1+2*K !; return
+           if (mod(tau+2,2)/=0) gamma = 2+2*K !; return
+           !
+         else
+           !
+           if (mod(K+tau+2,2)==0) gamma = 1+2*K  !; return
+           if (mod(K+tau+2,2)/=0) gamma = 2+2*K  !; return
+           !
+         endif 
+         !
       case('DNH','DNH(M)')
          !
          gamma = 0 

molecules.f90

@@ -818,8 +818,9 @@ recursive subroutine MLextF_potential(rank,ncoords,natoms,local,xyz,mu)
    real(ark)              :: force(molec%parmax) ! force(rank)
      !
      do i = 1,rank
-       !
-       if (extF%ifit(1,i)==0) then 
+       ! 
+       ! negative means the value of the parameter is fixed to the input value
+       if (extF%ifit(1,i)<0) then 
          force(i) = extF%coef(1,i)
        else
          force(i) = 0
@@ -834,8 +835,9 @@ recursive subroutine MLextF_potential(rank,ncoords,natoms,local,xyz,mu)
      do imu = 1,rank
        !
        !force = 0 
+       ! parameter with zero or negative indices are skipped from derivatives and expansions
        !
-       if (extF%ifit(1,imu)==0) cycle
+       if (extF%ifit(1,imu)<=0) cycle
        !
        force(imu) = 1.0_ark
        !
@@ -844,7 +846,7 @@ recursive subroutine MLextF_potential(rank,ncoords,natoms,local,xyz,mu)
        force(imu) = 0.0_ark
        !
      enddo
-    !
+     !
   end subroutine MLextF_potential
 
 

perturbation.f90

@@ -3575,6 +3575,12 @@ subroutine PTcontracted_prediagonalization(j)
                  cf%gamma = trim(sym%label(cf%isym))
                  !
                endif
+               if (trim(trove%symmetry)=='CSN') then
+                 !
+                 cf%isym = gamma + 2*contr(iclasses)%eigen(ilevel)%lquant
+                 cf%gamma = trim(sym%label(cf%isym))
+                 !
+               endif
                !
              endif
              !
@@ -32567,7 +32573,7 @@ subroutine PThamiltonianMat(jrot,nroots,diagonalizer_,uv_syevr_,postprocess_)
        write(out,"(/'Primitive matrix elements calculations...')")
     endif
     !
-    !$omp parallel do private(i,j,nu_i,nu_j,mat_elem) shared(a) schedule(static)
+    !$omp parallel do private(i,j,nu_i,nu_j,v_i,k_i,n_i,v_j,k_j,n_j,mat_elem) shared(a,b) schedule(dynamic)
     do i = 1,dimen
       !
       if (job%verbose>=5.and.mod(i,100)==0) print("('  i = ',i8)"), i