new O3 PES, modification of Legendre for negative sin(rho)), new fram…

…e for CH3OH
Trovemaster · Mar 1, 2024 · 68ae3a3 · 68ae3a3
1 parent 9ec4ef8
commit 68ae3a3
Show file tree

Hide file tree

Showing 5 changed files with 164 additions and 21 deletions.
diff --git 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
@@ -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 = <vl|h|vr> matrix elements and check if Numerov cracked
+              ! Now we test the h_t = <vl|h|vr> 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
@@ -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<small_) t2 = t2 + 2.0_ark*pi
+        if (t3-t2<small_) t3 = t3 + 2.0_ark*pi
+        !
+        theta12 = mod(t2-t1+2.0_ark*pi,2.0_ark*pi)
+        theta23 = mod(t3-t2+2.0_ark*pi,2.0_ark*pi)
+        theta13 = mod(t1-t3+2.0_ark*pi,2.0_ark*pi)
+        !
+        a1  = ( 2.0_ark*theta23 - theta13 - theta12 )/sqrt(6.0_ark)
+        a2  = (                   theta13 - theta12 )/sqrt(2.0_ark)
+        !
+        tbar = (t1 + t2 + t3-2.0_ark*pi)/3.0_ark
+        !
+        dst(10) = a1
+        dst(11) = a2
+        dst(12) = tbar 
+        !
+      else !  transform from TROVE coords to Z-matrix coords
+        !
+        dst(1:12) = src(1:12)+molec%local_eq(1:12)
+        !
+        A1 = src(10) 
+        A2 = src(11) 
+        tbar = src(12) ! + 2.0_ark*pi/3.0_ark
+        !
+        t1 = tbar+1.0_ark/3.0_ark*sqrt(2.0_ark)*A2 
+        t2 = 2.0_ark/3.0_ark*Pi+tbar-1.0_ark/6.0_ark*sqrt(2.0_ark)*A2-1.0_ark/6.0_ark*sqrt(6.0_ark)*A1
+        t3 = 4.0_ark/3.0_ark*Pi+tbar+1.0_ark/6.0_ark*sqrt(6.0_ark)*A1-1.0_ark/6.0_ark*sqrt(2.0_ark)*A2
+        !
+        dst(10) =  mod(t1+4.0_ark*pi,4.0_ark*pi)
+        dst(11) =  mod(t2+4.0_ark*pi,4.0_ark*pi)
+        dst(12) =  mod(t3+4.0_ark*pi,4.0_ark*pi)
+        !
+      endif
+      !
     end select
     !
     !
@@ -1492,9 +1540,9 @@ subroutine ML_b0_ch3oh(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
       tau = molec%taueq(1)
       !
       call ML_CH3OH_MEP_geometry(tau,r_eq(1),r_eq(2),r_eq(6),r_eq(3),r_eq(4),r_eq(5),&
-                                                beta1,beta2,beta3,alpha1,alpha2,alpha3,tau1,tau2,tau3,chi1,chi2,chi3)
+                                                 beta1,beta2,beta3,alpha1,alpha2,alpha3,tau1,tau2,tau3,chi1,chi2,chi3)
       !
-     if (trim(molec%potentype)=='POTEN_CH3OH_REF') then 
+      if (trim(molec%potentype)=='POTEN_CH3OH_REF') then 
 
 
         select case(trim(molec%coords_transform))

diff --git a/molecules.f90 b/molecules.f90
@@ -348,6 +348,10 @@ subroutine MLdefine_potenfunc
          !
          MLpotentialfunc => 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
@@ -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