diff --git a/fields.f90 b/fields.f90 index c9779d2..0b71812 100644 --- a/fields.f90 +++ b/fields.f90 @@ -18128,16 +18128,6 @@ subroutine FLbset1DNew(ibs,BSsize) stop 'At ME_laguerre_k: illegal imin' endif ! - if (f_t<=sqrt(small_).or.f_t>10000.0_rk) then - ! - write(out,"('1D ME_laguerre_k: the f_m constant is negative or too large.')") - write(out,"('Check if 1D PEC is quadratic at linearity if it is a linear molecule.')") - write(out,"('Alternatively, set f_m as a SPECPARAM parameter using sensible values from successful runs:')") - write(out,"('fm 0 77.0')") - stop '1D ME_laguerre_k: the f_m constant is negative or too large' - ! - endif - ! if (trove%specparam(bs%mode(1))>0d0 ) then ! f_t = trove%specparam(bs%mode(1))**2*g_t @@ -18149,6 +18139,16 @@ subroutine FLbset1DNew(ibs,BSsize) ! endif ! + if (f_t<=sqrt(small_).or.f_t>10000.0_rk) then + ! + write(out,"('1D ME_laguerre_k: the f_m constant is negative or too large.')") + write(out,"('Check if 1D PEC is quadratic at linearity if it is a linear molecule.')") + write(out,"('Alternatively, set f_m as a SPECPARAM parameter 3 using sensible values from successful runs:')") + write(out,"('fm 0 77.0')") + stop '1D ME_laguerre_k: the f_m constant is negative or too large' + ! + endif + ! f_m = sqrt(f_t/g_t) ! select case (trim(bs%type)) diff --git a/me_bnd.f90 b/me_bnd.f90 index 0b2b7f7..58cdafe 100644 --- a/me_bnd.f90 +++ b/me_bnd.f90 @@ -2288,7 +2288,7 @@ subroutine ME_Associate_Legendre(vmax,kmax,maxorder,rho_b_,isingular,npoints,drh ! rho = rho_b(1)+real(i,kind=ark)*rhostep x(i) = cos(rho) - sinrho(i) = sin(rho) + sinrho(i) = abs(sin(rho)) cosrho(i) = cos(rho) ! enddo @@ -2372,7 +2372,7 @@ subroutine ME_Associate_Legendre(vmax,kmax,maxorder,rho_b_,isingular,npoints,drh do i=0,npoints ! rho = rho_b(1)+real(i,kind=ark)*rhostep - phil(i) = L(i,vl)*sqrt(sin(rho)) + phil(i) = L(i,vl)*sqrt(abs(sin(rho))) dphil(i) = dL(i,vl)*sqrt(real((vl+k)*(vl-k+1),ark)) ! enddo @@ -2384,7 +2384,7 @@ subroutine ME_Associate_Legendre(vmax,kmax,maxorder,rho_b_,isingular,npoints,drh do i=0,npoints ! rho = rho_b(1)+real(i,kind=ark)*rhostep - phir(i) = L(i,vr)*sqrt(sin(rho)) + phir(i) = L(i,vr)*sqrt(abs(sin(rho))) dphir(i) = dL(i,vr)*sqrt(real((vr+k)*(vr-k+1),ark)) ! enddo @@ -2571,7 +2571,7 @@ subroutine ME_Associate_Legendre(vmax,kmax,maxorder,rho_b_,isingular,npoints,drh sigma = max(sigma,sigma_t) rms = rms + sigma_t**2 ! - ! Now we test the h_t = matrix elements and check if Numerov cracked + ! Now we test the h_t = matrix elements and check it cracked ! the Schroedinger all right if (nl/=nr.and.abs(h_t)>sqrt(small_)*abs(characvalue)*1e4) then write(out,"('ME_Associate_Legendre: wrong solution for <',i4,'|H|',i4,'> = ',f20.10)") nl,nr,h_t diff --git a/mol_ch3oh.f90 b/mol_ch3oh.f90 index 3800f32..f45e740 100644 --- a/mol_ch3oh.f90 +++ b/mol_ch3oh.f90 @@ -32,6 +32,7 @@ function ML_coordinate_transform_ch3oh(src,ndst,direct) result (dst) integer(ik) :: i0,nsrc, k(0:4),n(0:4),icoord real(ark) :: beta(3),alpha(3),cosa(3),phi1,phi2,phi3,chi1,chi2,chi3 real(ark) :: alpha0(3),cosa0(3),phi(3),chi(3),coschi(3) + real(ark) :: t1,t2,t3,theta12,theta23,theta13,a1,a2,tbar ! if (verbose>=7) write(out,"('ML_coordinate_transform_ch3oh/start')") ! @@ -679,8 +680,7 @@ function ML_coordinate_transform_ch3oh(src,ndst,direct) result (dst) !dst(10:12) = mod(dst(10:12)+2.0_ark*pi,2.0_ark*pi)+molec%local_eq(10:12) ! endif - - + ! case('R-ALPHA-S-TAU-BETA-REF') ! rref( 1,0:4) = molec%force( 1: 5) @@ -866,6 +866,54 @@ function ML_coordinate_transform_ch3oh(src,ndst,direct) result (dst) ! endif ! + case('R-ALPHA-THETA-TAU') + ! + if (direct) then + ! + !for stretches and 'alpha' bends just subtract equilibrium coordinates + dst(1:12) = src(1:12)-molec%local_eq(1:12) + ! + t1 = src(10) + t2 = src(11) + t3 = src(12) + ! + ! subtract equilbrium theta values to make a1/a2 zero at equilibrium + ! and ensure consistent transfroms + ! + if (t2-t1 MLpoten_c3_R_theta ! + case('POTEN_O3_OLEG') + ! + MLpotentialfunc => MLpoten_o3_oleg + ! case('POTEN_H2S_DVR3D') ! MLpotentialfunc => MLpoten_h2s_dvr3d @@ -490,7 +494,7 @@ subroutine MLdefine_kinetic_subroutine ! MLkineticfunc => MLkinetic_xyz_bond_EKE ! - case('KINETIC_XYZ_EKE_BISECT_SINRHO') + case('KINETIC_XYZ_EKE_BOND_SINRHO') ! MLkineticfunc => MLkinetic_xyz_EKE_sinrho ! diff --git a/pot_xy2.f90 b/pot_xy2.f90 index e2412a4..2c38505 100644 --- a/pot_xy2.f90 +++ b/pot_xy2.f90 @@ -4,7 +4,6 @@ module pot_xy2 use accuracy use moltype - !use pot_h216o implicit none @@ -16,7 +15,8 @@ module pot_xy2 MLdipole_xy2_lorenzo,MLdms2pqr_xy2_linear,MLpoten_xy2_bubukina,MLpoten_xyz_tyuterev,MLdms2pqr_xyz_coeff,& MLpoten_xy2_tyuterev_alpha,MLdms2pqr_xyz_z_frame,MLpoten_xyz_tyuterev_rho public MLpoten_xy2_morse_cos,MLpoten_xyz_Koput,MLdipole_bisect_s1s2theta_xy2,MLloc2pqr_xyz,MLpoten_xy2_sym_morse - public MLdms_c3_Schroeder,MLpoten_xy2_morse_powers,MLdms2pqr_xyz_z_frame_sinrho,MLdms_Schroeder_xyz_Eckart + public MLdms_c3_Schroeder,MLpoten_xy2_morse_powers,MLdms2pqr_xyz_z_frame_sinrho,MLdms_Schroeder_xyz_Eckart,& + MLpoten_o3_oleg private integer(ik), parameter :: verbose = 4 ! Verbosity level @@ -6972,6 +6972,97 @@ function MLpoten_xy2_morse_powers(ncoords,natoms,local,xyz,force) result(f) end function MLpoten_xy2_morse_powers + ! Ozone potential function from Polyansky et al 2017 + ! JQSRT doi:10.1016/j.jqsrt.2018.02.018 + ! + function MLpoten_o3_oleg(ncoords,natoms,local,xyz,force) result(f) + ! + integer(ik),intent(in) :: ncoords,natoms + real(ark),intent(in) :: local(ncoords) + real(ark),intent(in) :: xyz(natoms,3) + real(ark),intent(in) :: force(:) + real(ark) :: f + ! + real(ark) :: r1,r2,theta,re,thetae,e0,a1,b1,b2,c1 + real(ark) :: xs1,xs2,xst,vpot,v0,v1,roo,voo,xep1,xep2,xep3,xep4,rr1,rr2 + real(ark) :: pd,deg + ! + if (verbose>=6) write(out,"('MLpoten_o3_oleg/start')") + ! + pd=pi/180.0_ark + deg=180.0_ark/pi + ! + r1 = local(1) + r2 = local(2) + theta = local(3) + ! + re = force(1) + thetae = force(2)*pd + e0 = force(3) + a1 = force(4) + b1 = force(5) + b2 = force(6) + c1 = force(7) + ! + xs1 = (r1+r2)*0.5_ark-re + xs2 = (r1-r2)*0.5_ark + xst = cos(theta)-cos(thetae) + ! + rr1 = r1-re + rr2 = r2-re + roo = sqrt(r1**2+r2**2-2.0_ark*r1*r2*cos(theta)) + voo = 800000.0_ark*exp(-c1*roo) + ! + xep1 = (exp(-2.0_ark*a1*rr1)-2.0_ark*exp(-a1*rr1)+1.0_ark)*e0 + xep2 = (exp(-2.0_ark*a1*rr2)-2.0_ark*exp(-a1*rr2)+1.0_ark)*e0 + xep3 = exp(-b1*((r1-re)**2+(r2-re)**2)) + xep4 = exp(-b2*(theta-thetae)) + ! + v0 = force( 8)**xs1**0*xs2**0*xst**0 + v1 = force( 9)*xs1**0*xs2**0*xst**1& + +force( 10)*xs1**1*xs2**0*xst**0& + +force( 11)*xs1**0*xs2**0*xst**2& + +force( 12)*xs1**0*xs2**2*xst**0& + +force( 13)*xs1**1*xs2**0*xst**1& + +force( 14)*xs1**2*xs2**0*xst**0& + +force( 15)*xs1**0*xs2**0*xst**3& + +force( 16)*xs1**0*xs2**2*xst**1& + +force( 17)*xs1**1*xs2**0*xst**2& + +force( 18)*xs1**1*xs2**2*xst**0& + +force( 19)*xs1**2*xs2**0*xst**1& + +force( 20)*xs1**3*xs2**0*xst**0& + +force( 21)*xs1**0*xs2**0*xst**4& + +force( 22)*xs1**0*xs2**2*xst**2& + +force( 23)*xs1**0*xs2**4*xst**0& + +force( 24)*xs1**1*xs2**0*xst**3& + +force( 25)*xs1**1*xs2**2*xst**1& + +force( 26)*xs1**2*xs2**0*xst**2& + +force( 27)*xs1**2*xs2**2*xst**0& + +force( 28)*xs1**3*xs2**0*xst**1& + +force( 29)*xs1**4*xs2**0*xst**0& + +force( 30)*xs1**0*xs2**0*xst**5& + +force( 31)*xs1**0*xs2**2*xst**3& + +force( 32)*xs1**0*xs2**4*xst**1& + +force( 33)*xs1**1*xs2**0*xst**4& + +force( 34)*xs1**1*xs2**2*xst**2& + +force( 35)*xs1**1*xs2**4*xst**0& + +force( 36)*xs1**2*xs2**0*xst**3& + +force( 37)*xs1**2*xs2**2*xst**1& + +force( 38)*xs1**3*xs2**0*xst**2& + +force( 39)*xs1**3*xs2**2*xst**0& + +force( 40)*xs1**4*xs2**0*xst**1& + +force( 41)*xs1**5*xs2**0*xst**0& + +force( 42)*xs1**0*xs2**0*xst**6& + +force( 43)*xs1**0*xs2**6*xst**0 + ! + vpot=v0+v1*xep3+voo+xep1+xep2+xep4 + if(theta*deg< 50.0_ark) vpot=30000.0_ark + ! + f = vpot + ! + if (verbose>=6) write(out,"('MLpoten_o3_oleg/end')") + ! + end function MLpoten_o3_oleg