From c9d15f6e5dab1eb991060446ad583ee28d03179f Mon Sep 17 00:00:00 2001 From: Sergey Yurchenko Date: Fri, 25 Aug 2023 12:33:52 +0100 Subject: [PATCH] CSNV, bug in contracted basis for sing_triatom, compact form for fields and energies, fix, fit keywords, H2CO PES with isotopologues --- fields.f90 | 148 +++++++++++++++++++++++++++++++++-------------- kin_xy2.f90 | 4 +- mol_xy2.f90 | 36 +++++++++++- molecules.f90 | 10 ++-- perturbation.f90 | 8 ++- pot_zxy2.f90 | 60 ++++++++++++++----- refinement.f90 | 35 ++++++++--- symmetry.f90 | 4 +- 8 files changed, 231 insertions(+), 74 deletions(-) diff --git a/fields.f90 b/fields.f90 index 3f65f9d..d7c44f3 100644 --- a/fields.f90 +++ b/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.nitemstrove%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) diff --git a/kin_xy2.f90 b/kin_xy2.f90 index a3b130e..ca6fbb6 100644 --- a/kin_xy2.f90 +++ b/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) diff --git a/mol_xy2.f90 b/mol_xy2.f90 index 10a381f..132afdb 100644 --- a/mol_xy2.f90 +++ b/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 diff --git a/molecules.f90 b/molecules.f90 index 188fac2..750f91a 100644 --- a/molecules.f90 +++ b/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 diff --git a/perturbation.f90 b/perturbation.f90 index 725fee8..84d4cb4 100644 --- a/perturbation.f90 +++ b/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 diff --git a/pot_zxy2.f90 b/pot_zxy2.f90 index c00e6b9..56d9643 100644 --- a/pot_zxy2.f90 +++ b/pot_zxy2.f90 @@ -1298,6 +1298,7 @@ function MLpoten_zxy2_mep_r_alpha_rho_powers(ncoords,natoms,x,xyz,force) result( end function MLpoten_zxy2_mep_r_alpha_rho_powers + function MLpoten_zxy2_mep_r_alpha_rho_powers_iso(ncoords,natoms,x,xyz,force) result(f) ! integer(ik),intent(in) :: ncoords,natoms @@ -1306,13 +1307,16 @@ function MLpoten_zxy2_mep_r_alpha_rho_powers_iso(ncoords,natoms,x,xyz,force) res real(ark),intent(in) :: force(:) real(ark) :: f ! - integer(ik) :: N,N0,N_iso,iterm,k_ind(6) - real(ark) :: M0,M_main,V_iso,V,M_iso,xieq(6),local(6),cosrho,rho,y(6),xi(6) + integer(ik) :: N,N0,N_iso_carbon,N_iso_hydrogen,iterm,k_ind(6) + real(ark) :: M0,CarbonMassMain,HydrogenMassMain,V_iso_carbon,V_iso_hydrogen,V,CarbonMassIso,HydrogenMassIso,xieq(6),local(6),cosrho,rho,y(6),xi(6) ! N0 = force(1) - N_iso = force(2) - M_main = force(3) - M_iso = force(4) + N_iso_carbon = force(2) + N_iso_hydrogen = force(3) + CarbonMassMain = force(4) + CarbonMassIso = force(5) + HydrogenMassMain = force(6) + HydrogenMassIso = force(7) ! local = from_local_to_r1r2r3a1a2tau(x,6) ! @@ -1324,13 +1328,13 @@ function MLpoten_zxy2_mep_r_alpha_rho_powers_iso(ncoords,natoms,x,xyz,force) res ! !xieq(1:6) = ML_MEP_zxy2_R_rho(rho) ! - xieq(1) = sum(force(5:9)*cosrho**molec%pot_ind(1,5:9)) - xieq(2) = sum(force(10:14)*cosrho**molec%pot_ind(2,10:14)) + xieq(1) = sum(force(8:12)*cosrho**molec%pot_ind(1,8:12)) + xieq(2) = sum(force(13:17)*cosrho**molec%pot_ind(2,13:17)) xieq(3) = xieq(2) - xieq(4) = sum(force(15:19)*cosrho**molec%pot_ind(4,15:19)) + xieq(4) = sum(force(18:22)*cosrho**molec%pot_ind(4,18:22)) xieq(5) = xieq(4) ! - N = 19 + N = 22 ! ! expansion functions ! @@ -1366,13 +1370,13 @@ function MLpoten_zxy2_mep_r_alpha_rho_powers_iso(ncoords,natoms,x,xyz,force) res enddo N = N+N0 ! - V_iso = 0.0_ark + V_iso_carbon = 0.0_ark ! - do iterm = N+1,N+N_iso + do iterm = N+1,N+N_iso_carbon ! xi(1:6) = y(1:6)**molec%pot_ind(1:6,iterm) ! - V_iso = V_iso + force(iterm)*product(xi(1:6)) + V_iso_carbon = V_iso_carbon + force(iterm)*product(xi(1:6)) ! if (molec%pot_ind(2,iterm)/=molec%pot_ind(3,iterm).or.molec%pot_ind(4,iterm)/=molec%pot_ind(5,iterm)) then ! @@ -1385,16 +1389,42 @@ function MLpoten_zxy2_mep_r_alpha_rho_powers_iso(ncoords,natoms,x,xyz,force) res ! xi(1:6) = y(1:6)**k_ind(1:6) ! - V_iso = V_iso + force(iterm)*product(xi(1:6)) + V_iso_carbon = V_iso_carbon + force(iterm)*product(xi(1:6)) ! endif ! enddo ! - f = V + V_iso*(M_main-M_iso)/M_main + N = N+N_iso_carbon ! -end function MLpoten_zxy2_mep_r_alpha_rho_powers_iso + V_iso_hydrogen = 0.0_ark + ! + do iterm = N+1,N+N_iso_hydrogen + ! + xi(1:6) = y(1:6)**molec%pot_ind(1:6,iterm) + ! + V_iso_hydrogen = V_iso_hydrogen + force(iterm)*product(xi(1:6)) + ! + if (molec%pot_ind(2,iterm)/=molec%pot_ind(3,iterm).or.molec%pot_ind(4,iterm)/=molec%pot_ind(5,iterm)) then + ! + k_ind(1) = molec%pot_ind(1,iterm) + k_ind(2) = molec%pot_ind(3,iterm) + k_ind(3) = molec%pot_ind(2,iterm) + k_ind(4) = molec%pot_ind(5,iterm) + k_ind(5) = molec%pot_ind(4,iterm) + k_ind(6) = molec%pot_ind(6,iterm) + ! + xi(1:6) = y(1:6)**k_ind(1:6) + ! + V_iso_hydrogen = V_iso_hydrogen + force(iterm)*product(xi(1:6)) + ! + endif + ! + enddo + f = V + V_iso_carbon*(CarbonMassIso-CarbonMassMain)/CarbonMassMain + V_iso_hydrogen*(HydrogenMassIso - HydrogenMassMain)/HydrogenMassMain + ! +end function MLpoten_zxy2_mep_r_alpha_rho_powers_iso ! ! Defining potential energy function diff --git a/refinement.f90 b/refinement.f90 index 8bdaa63..a5f15b5 100644 --- a/refinement.f90 +++ b/refinement.f90 @@ -371,7 +371,7 @@ subroutine sf_fitting(Jval) character(len=1) :: rng character(len=1),allocatable :: mark(:) character(len=cl) :: my_fmt,my_fmt_pot1,my_fmt_pot2 !format for I/O specification - character(len=cl) :: my_fmt_en2,my_fmt_par1,my_fmt_par2 !format for I/O specification + character(len=cl) :: my_fmt_en2,my_fmt_par1,my_fmt_par2,my_fmt_par_fit !format for I/O specification character(len=wl) :: my_fmt_en1 !wider format for I/O specification ! if (job%verbose>=2) write(out,"(/'The least-squares fitting ...')") @@ -592,7 +592,8 @@ subroutine sf_fitting(Jval) write(my_fmt_pot1,'(a,i0,a)') "(1h1,5x,a,a,a//4x,",ncoords,"(7x),a,7x,a,3x,a,3x,a/)" write(my_fmt_pot2,'(a1,i0,a)') "(",ncoords,"(2x,f18.9),2(1x,g12.5),1x,f12.5,1x,e12.4)" - write(my_fmt_par1,'(a,i0,a)') "(a8,4x,",Ncoords,"i3,1x,i2,e22.14)" + write(my_fmt_par1, '(a,i0,a)') "(a8,4x,",Ncoords,"i3,2x,e22.14)" + write(my_fmt_par_fit,'(a,i0,a)') "(a8,4x,",Ncoords,"i3,2x,e22.14,4x,a3)" ! nlevels = Neigenlevels ! @@ -1397,7 +1398,7 @@ subroutine sf_fitting(Jval) ! do ncol=1,numpar i = ifitparam(ncol) - potparam(i)=potparam(i)+dx(ncol)*fitting%fit_scale + if (ivar(i)>0) potparam(i)=potparam(i)+dx(ncol)*fitting%fit_scale enddo ! ! Robust fit: adjust the fitting weights @@ -1468,13 +1469,23 @@ subroutine sf_fitting(Jval) if (all(extF%term(:,1,1)==-1)) then ! do i=1,parmax - write (out,"(a8,4x,i2,e22.14)") nampar(i),ivar(i),potparam(i) + if (ivar(i)>0) then + write (out,"(a8,4x,e22.14,4x,'fit')") nampar(i),potparam(i) + elseif (ivar(i)<0) then + write (out,"(a8,4x,e22.14,4x,'fix')") nampar(i),potparam(i) + else + write (out,"(a8,4x,e22.14)") nampar(i),potparam(i) + endif enddo ! write(out,"(/'Potential parameters:')") ! do i=1,parmax - write (out,"(a8,4x,i2,e22.14)") nampar(i),ivar(i),molec%force(i)+potparam(i) + if (ivar(i)>0) then + write (out,"(a8,4x,e22.14)") nampar(i),molec%force(i)+potparam(i) + else + write (out,"(a8,4x,e22.14)") nampar(i),molec%force(i) + endif enddo ! else @@ -1482,13 +1493,23 @@ subroutine sf_fitting(Jval) ! 'powers' ! do i=1,parmax - write (out,my_fmt_par1) nampar(i),(extF%term(l,1,i),l=1,Ncoords),ivar(i),potparam(i) + if (ivar(i)>0) then + write (out,my_fmt_par_fit) nampar(i),(extF%term(l,1,i),l=1,Ncoords),potparam(i),"fit" + elseif (ivar(i)<0) then + write (out,my_fmt_par_fit) nampar(i),(extF%term(l,1,i),l=1,Ncoords),potparam(i),"fix" + else + write (out,my_fmt_par1) nampar(i),(extF%term(l,1,i),l=1,Ncoords),potparam(i) + endif enddo ! write(out,"(/'Potential parameters:')") ! do i=1,parmax - write (out,my_fmt_par1) nampar(i),(extF%term(l,1,i),l=1,Ncoords),ivar(i),molec%force(i)+potparam(i) + if (ivar(i)>0) then + write (out,my_fmt_par1) nampar(i),(extF%term(l,1,i),l=1,Ncoords),molec%force(i)+potparam(i) + else + write (out,my_fmt_par1) nampar(i),(extF%term(l,1,i),l=1,Ncoords),molec%force(i) + endif enddo ! endif diff --git a/symmetry.f90 b/symmetry.f90 index a2cc745..9e9b0a9 100644 --- a/symmetry.f90 +++ b/symmetry.f90 @@ -264,7 +264,7 @@ subroutine SymmetryInitialize(sym_group) elseif(mod(j,2) == 0) then sym_sub_label = 'f' endif - write(k_num,'(i5)') (j-1-mod(j-1,4))/4 + write(k_num,'(i5)') (j-1-mod(j-1,2))/2 sym_cur_label = trim(sym_sub_label)//trim(adjustl(k_num)) sym%label(j) = sym_cur_label !write(*,*) sym%label(j) @@ -280,7 +280,7 @@ subroutine SymmetryInitialize(sym_group) sym%product_table(1:2,1:2) = reshape((/0,0, & 0,0/), (/2,2/)) ! - do j = 1,2**(n_cs-2)-1 + do j = 1,2**(n_cs-1)-1 sym%product_table(2*j+1: 2*j+2,1:2) = reshape((/1,2, & 1,1 /), (/2,2/)) enddo