Skip to content

Commit

Permalink
Lots of changes
Browse files Browse the repository at this point in the history 
MLpoten_xyz_morse_fourier_mep
CH3OH 'R-ALPHA-THETA-TAU'
ZXY3 print b0 cartesian
CH4 module adding CH3D r-beta-sym and  C3V
me_bnd laguerre realtive ZPE
fields laguerre-k redefining fm for f_t larger
  • Loading branch information
@Trovemaster
Trovemaster committed May 15, 2024
1 parent ea1953d commit b2a7ef8
Show file tree
Hide file tree
Showing 10 changed files with 1,001 additions and 309 deletions.
2 changes: 1 addition & 1 deletion docs/source/quantumnumbers.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -301,7 +301,7 @@ The ro-vibrational calculations in the :math:`J=0` representation (step 3) are c
Phi_{\lambda}^{(J=0)} |J,K,\tau\rangle,
where :math`| J,K,\Gamma_{\rm rot} \rangle ` is a rigid rotor basis function. Therefore, in this case the energy output has a slightly different form:
where :math:`| J,K,\Gamma_{\rm rot} \rangle` is a rigid rotor basis function. Therefore, in this case the energy output has a slightly different form:
::

------ ------ ------------- ------- --- -- --- ----- ----- --- ---- -------- ------ --- --- --- --------
Expand Down
4 changes: 2 additions & 2 deletions fields.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18099,7 +18099,7 @@ subroutine FLbset1DNew(ibs,BSsize)
!
! estimate second derivative using a step of ~0.005 rad
!
! irho that that distance:
! irho at that distance:
irho = max(nint(fd_step/trove%rhostep),2)
!
f_t = ( 2.0_ark*f1drho(irho)-2.0_ark*f1drho(0) )/(fd_step)**2
Expand Down Expand Up @@ -18139,7 +18139,7 @@ subroutine FLbset1DNew(ibs,BSsize)
!
endif
!
if (f_t<=sqrt(small_).or.f_t>10000.0_rk) then
if (f_t<=sqrt(small_).or.f_t>500000.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.')")
Expand Down
4 changes: 2 additions & 2 deletions me_bnd.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2216,7 +2216,7 @@ subroutine ME_Associate_Legendre(vmax,kmax,maxorder,rho_b_,isingular,npoints,drh
integer(ik),intent(in) :: icoord ! coordinate number for which the numerov is employed
integer(ik),intent(in) :: verbose ! Verbosity level
!
real(ark) :: rho,rhostep,potmin,C_l,C_r,ddphi_rho(vmax+1),zpe
real(ark) :: rho,rhostep,potmin,C_l,C_r,ddphi_rho(vmax+1),zpe = 0
real(ark) :: psipsi_t,characvalue,rho_b(2),h_t,sigma_t,sigma,rms,C1,C2,C3,C4,cross_prod,factor,mu_zz_t,mu_rr_t
!
integer(ik) :: vl,vr,nl,nr,il,ir,nmax,lambda,alloc,i,k,rec_len,n,imin,io_slot,lmax,nmax1
Expand Down Expand Up @@ -2439,7 +2439,7 @@ subroutine ME_Associate_Legendre(vmax,kmax,maxorder,rho_b_,isingular,npoints,drh
!
write (out,"(/' Legendre-optimized energies are:')")
!
zpe = ener(1)
if (k==0) zpe = ener(1)
!
do nl=0,nmax
i = k*(nmax+1)+nl
Expand Down
188 changes: 90 additions & 98 deletions mol_ch3oh.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1514,51 +1514,44 @@ subroutine ML_b0_ch3oh(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
!
tau = molec%taueq(1)
!
if (trim(molec%potentype)=='POTEN_CH3OH_REF') then


!
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)
!
select case(trim(molec%coords_transform))
case default
write (out,"('ML_b0_ch3oh: coord. type ',a,' unknown')") trim(molec%coords_transform)
stop 'ML_b0_ch3oh - bad coord. type'
!
case('R-TAU-BOWMAN-REF')

!if (trim(molec%potentype)=='POTEN_CH3OH_REF'.or.trim(molec%potentype)=='R-ALPHA-S-TAU-BETA-REF'&
!.or.trim(molec%potentype)=='R-TAU-BOWMAN-REF'.or.trim(molec%potentype)=='R-1TAU-2BETA-REF') then

!
r_at(1:6) = r_eq(1:6)
!
r_at(7) = beta1
r_at(8) = beta2
r_at(9) = alpha2
r_at(10) = beta3
r_at(11) = alpha3
r_at(12) = tau
!
case('R-TAU-CHI')
!
r_at(1:6) = r_eq(1:6)
!
r_at(7) = beta1
r_at(8) = beta2
r_at(9) = beta3
r_at(10) = tau1
r_at(11) = tau2
r_at(12) = tau3
!
end select
!
else
!
r_at = molec%local_eq
!
endif
r_at = molec%local_eq
!
!
!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)
!
select case(trim(molec%frame))
case default
write (out,"('ML_b0_ch3oh: coord. type ',a,' unknown')") trim(molec%coords_transform)
stop 'ML_b0_ch3oh - bad coord. type'
!
case('R-TAU-BOWMAN-REF')
!
r_at(1:6) = r_eq(1:6)
!
r_at(7) = beta1
r_at(8) = beta2
r_at(9) = alpha2
r_at(10) = beta3
r_at(11) = alpha3
r_at(12) = tau
!
case('R-TAU-CHI')
!
r_at(1:6) = r_eq(1:6)
!
r_at(7) = beta1
r_at(8) = beta2
r_at(9) = beta3
r_at(10) = tau1
r_at(11) = tau2
r_at(12) = tau3
!
case('R-ALPHA-THETA-TAU')
!
r_at = molec%local_eq
!
end select
!
call MLfromlocal2cartesian(-1,r_at,a0_ark)
!
Expand All @@ -1585,60 +1578,57 @@ subroutine ML_b0_ch3oh(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
!
tau = rho_i(i)
!
if (trim(molec%potentype)=='POTEN_CH3OH_REF') then
!
!
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)
!
select case(trim(molec%coords_transform))
case default
write (out,"('ML_b0_ch3oh: coord. type ',a,' unknown')") trim(molec%coords_transform)
stop 'ML_b0_ch3oh - bad coord. type'
!
case('R-TAU-BOWMAN-REF')

!if (trim(molec%potentype)=='POTEN_CH3OH_REF'.or.trim(molec%potentype)=='R-ALPHA-S-TAU-BETA-REF'&
!.or.trim(molec%potentype)=='R-TAU-BOWMAN-REF'.or.trim(molec%potentype)=='R-1TAU-2BETA-REF') then

!
r_at(1:6) = r_eq(1:6)
!
r_at(7) = beta1
r_at(8) = beta2
r_at(9) = alpha2
r_at(10) = beta3
r_at(11) = alpha3
r_at(12) = tau
!
case('R-TAU-CHI')
!
r_at(1:6) = r_eq(1:6)
!
r_at(7) = beta1
r_at(8) = beta2
r_at(9) = beta3
r_at(10) = tau1
r_at(11) = tau2
r_at(12) = tau3
!
end select
!
!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)
!
select case(trim(molec%coords_transform))
case default
write (out,"('ML_b0_ch3oh: coord. type ',a,' unknown')") trim(molec%coords_transform)
stop 'ML_b0_ch3oh - bad coord. type'
!
case('R-TAU-BOWMAN-REF')

!if (trim(molec%potentype)=='POTEN_CH3OH_REF'.or.trim(molec%potentype)=='R-ALPHA-S-TAU-BETA-REF'&
!.or.trim(molec%potentype)=='R-TAU-BOWMAN-REF'.or.trim(molec%potentype)=='R-1TAU-2BETA-REF') then

!
r_at(1:6) = r_eq(1:6)
!
r_at(7) = beta1
r_at(8) = beta2
r_at(9) = alpha2
r_at(10) = beta3
r_at(11) = alpha3
r_at(12) = tau
!
case('R-TAU-CHI')
!
if (verbose>=6) then
write(out,"('r0',i8,12f12.8)") i,r_at(1:12)
endif
r_at(1:6) = r_eq(1:6)
!
else
r_at(7) = beta1
r_at(8) = beta2
r_at(9) = beta3
r_at(10) = tau1
r_at(11) = tau2
r_at(12) = tau3
!
case('R-ALPHA-THETA-TAU')
!
r_eq(10) = tau
r_eq(11) = tau+molec%local_eq(10)-molec%local_eq(11)
r_eq(12) = tau+molec%local_eq(10)-molec%local_eq(12)
!
r_at = r_eq
!
end select
!
if (verbose>=6) then
write(out,"('r0',i8,12f12.8)") i,r_at(1:12)
endif
!
call MLfromlocal2cartesian(-1,r_at,a0_ark)
!
b0(:,:,i) = a0_ark(:,:)
!
call MLorienting_a0(molec%Natoms,molec%AtomMasses,b0(:,:,i),transform)
Expand All @@ -1655,25 +1645,27 @@ subroutine ML_b0_ch3oh(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
! write(out,"('b0',i8,18f12.8)") i,b0(:,:,i)
!endif

if (verbose>=5) then
enddo
!
if (verbose>=4) then
do i = 0,npoints
write(out,"(i6)") molec%natoms
!
write(out,"(/'C',3x,3f14.8)") b0(1,:,i)
write(out,"( 'O',3x,3f14.8)") b0(2,:,i)
write(out,"( 'H',3x,3f14.8)") b0(3,:,i)
write(out,"( 'H',3x,3f14.8)") b0(4,:,i)
write(out,"( 'H',3x,3f14.8)") b0(5,:,i)
write(out,"( 'H',3x,3f14.8)") b0(6,:,i)
endif

enddo
write(out,"(/a,3x,3f14.8)") trim(molec%zmatrix(1)%name),b0(1,:,i)
write(out,"( a,3x,3f14.8)") trim(molec%zmatrix(2)%name),b0(2,:,i)
write(out,"( a,3x,3f14.8)") trim(molec%zmatrix(3)%name),b0(3,:,i)
write(out,"( a,3x,3f14.8)") trim(molec%zmatrix(3)%name),b0(4,:,i)
write(out,"( a,3x,3f14.8)") trim(molec%zmatrix(3)%name),b0(5,:,i)
write(out,"( a,3x,3f14.8)") trim(molec%zmatrix(3)%name),b0(6,:,i)
!
enddo
endif
!
do i = 0,npoints
if (verbose>=6) then
write(out,"('b0',i8,18f12.8)") i,( (b0(i0,icoord,i),icoord=1,3),i0=1,6 )
endif
enddo

!
endif
!
Expand Down
Loading

0 comments on commit b2a7ef8

Please sign in to comment.