173 lines
4.8 KiB
Fortran
173 lines
4.8 KiB
Fortran
module pesmod
|
|
!use dim_parameter,only:qn,ndiab,pst
|
|
use iso_fortran_env, only: idp => int32, dp => real64
|
|
use dip_param
|
|
|
|
use nn_params
|
|
use nncommons, only: len_in
|
|
implicit none
|
|
|
|
integer(idp),parameter:: ndiab=4
|
|
private
|
|
public :: potential_nh3_low
|
|
contains
|
|
|
|
|
|
subroutine potential_nh3_low(e,q,nn_out)
|
|
implicit none
|
|
real(dp),intent(in)::q(maxnin)
|
|
real(dp),intent(out)::e(:,:)
|
|
real(dp),intent(inout):: nn_out(maxnout)
|
|
!logical,intent(in):: sym
|
|
!integer(idp),intent(in)::key
|
|
integer(idp) id,i, ii
|
|
real(dp) tmp,xs,xb,ys,yb,a,b,ss,sb
|
|
logical, parameter:: switch = .false. ! parameter to switch when you want to know Ex and Ey
|
|
logical, parameter:: sym = .false. ! to include r**2 in the Genetic fitting
|
|
real(dp),parameter:: tol = 1.0d-9
|
|
real(dp),dimension(16):: scal, shift
|
|
!real(dp),dimension(17):: shi,scalee
|
|
!include "shift-scale-70N.incl"
|
|
!include "sh-scal-50.incl"
|
|
xs=q(len_in+2)
|
|
ys=q(len_in+3)
|
|
xb=q(len_in+4)
|
|
yb=q(len_in+5)
|
|
a=q(len_in+1)
|
|
b=q(len_in+6)
|
|
|
|
ss= xs**2+ys**2
|
|
sb= xb**2+yb**2
|
|
call init_dip_planar_data()
|
|
shift =1.0d0
|
|
scal = 1.0d-3
|
|
!scal = 1.0d-3
|
|
!scal(1:3) = 1.0d-3
|
|
ii = 1
|
|
do i =1 ,27 ! pes parameter key index
|
|
if( abs(p(i)) .gt. tol ) then
|
|
p(i) = p(i)*(shift(ii)+ scal(ii)*nn_out(ii))
|
|
ii = ii + 1
|
|
|
|
endif
|
|
enddo
|
|
!write(6,*)" PES non zero", ii-1
|
|
|
|
|
|
e = 0.0_dp
|
|
id = 1 ! 1
|
|
! v -term
|
|
! oth order only
|
|
! A2"
|
|
e(1,1) = e(1,1) + p(pst(1,id)) !- p(pst(1,1))
|
|
if (sym) then
|
|
e(1,1)=e(1,1) + p(pst(1,id)+1)*ss +p(pst(1,id)+2)*sb
|
|
endif
|
|
id = id +1 !2
|
|
! E block
|
|
e(2,2) = e(2,2)+ p(pst(1,id))! - p(pst(1,1))
|
|
e(3,3) = e(3,3)+ p(pst(1,id))! - p(pst(1,1))
|
|
if (sym) then
|
|
e(2,2)=e(2,2)+p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2
|
|
e(3,3)=e(3,3)+p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2
|
|
endif
|
|
|
|
id = id +1 ! 3
|
|
|
|
!A1'
|
|
e(4,4) = e(4,4) + p(pst(1,id))! - p(pst(1,1))
|
|
if (sym) then
|
|
e(4,4)=e(4,4) + p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2 &
|
|
+ p(pst(1,id)+3)*ss**3
|
|
endif
|
|
|
|
! W and z of JT block
|
|
|
|
! first order
|
|
id = id +1 !4
|
|
|
|
e(2,2) = e(2,2) + p(pst(1,id))*xs + p(pst(1,id)+1)*xb
|
|
e(3,3) = e(3,3) - p(pst(1,id))*xs - p(pst(1,id)+1)*xb
|
|
e(2,3) = e(2,3) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb
|
|
|
|
! order 2
|
|
id = id +1 ! 5
|
|
e(2,2) = e(2,2) + p(pst(1,id))*(xs**2-ys**2) &
|
|
+ p(pst(1,id)+1)*(xb**2-yb**2) &
|
|
+ p(pst(1,id)+2)*(xs*xb - ys*yb)
|
|
e(3,3) = e(3,3) - p(pst(1,id))*(xs**2-ys**2) &
|
|
- p(pst(1,id)+1)*(xb**2-yb**2) &
|
|
- p(pst(1,id)+2)*(xs*xb - ys*yb)
|
|
e(2,3) = e(2,3) - p(pst(1,id))*(2*xs*ys) &
|
|
- p(pst(1,id)+1)*(2*xb*yb) &
|
|
- p(pst(1,id)+2)*(xs*yb+xb*ys)
|
|
|
|
! the coupling of A2' and E'
|
|
! order 1
|
|
id = id +1 !6
|
|
|
|
e(1,2) = e(1,2) + b*(p(pst(1,id))*xs + p(pst(1,id)+1)*yb)
|
|
e(1,3) = e(1,3) - b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
|
|
! order 2
|
|
id =id +1 ! 7
|
|
e(1,2) = e(1,2) +b*(p(pst(1,id))*(xs**2-ys**2) &
|
|
+ p(pst(1,id)+1)*(xb**2-yb**2) &
|
|
+ p(pst(1,id)+2)*(xs*xb -ys*yb))
|
|
e(1,3) = e(1,3) +b*(p(pst(1,id))*(2*xs*ys) &
|
|
+ p(pst(1,id)+1)*(2*xb*yb) &
|
|
+ p(pst(1,id)+2)*(xs*yb + xb*ys))
|
|
|
|
! the coupling of A1' and E'
|
|
id = id +1 ! 8
|
|
e(2,4) = e(2,4) + p(pst(1,id))*xs + p(pst(1,id)+1)*xb
|
|
e(3,4) = e(3,4) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb
|
|
if (sym) then
|
|
e(2,4) = e(2,4)+ p(pst(1,id)+2)*xs*ss
|
|
e(3,4) = e(3,4)- p(pst(1,id)+2)*ys*ss
|
|
endif
|
|
|
|
! order 2
|
|
id = id +1 ! 9
|
|
e(2,4) = e(2,4) + p(pst(1,id))*(xs**2-ys**2) &
|
|
+ p(pst(1,id)+1)*(xb**2-yb**2) &
|
|
+ p(pst(1,id)+2)*(xs*xb -ys*yb)
|
|
e(3,4) = e(3,4) + p(pst(1,id))*(2*xs*ys) &
|
|
+ p(pst(1,id)+1)*(2*xb*yb) &
|
|
+ p(pst(1,id)+2)*(xs*yb + xb*ys)
|
|
! the couplign between A
|
|
id = id +1 ! 10
|
|
e(1,4) = e(1,4) + b*p(pst(1,id))
|
|
!write(*,*) 'Id = ',id
|
|
|
|
call copy_2_lower_triangle(e)
|
|
if (switch) then
|
|
write( 6,*)" EX AND EY ARE SWITCHED"
|
|
tmp = e(2,2)
|
|
e(2,2) = e(3,3)
|
|
e(3,3) = tmp
|
|
|
|
! the coupling
|
|
|
|
tmp = e(2,4)
|
|
e(2,4) = e(3,4)
|
|
e(3,4) = tmp
|
|
e(4,2) = e(2,4)
|
|
e(4,3) = e(3,4)
|
|
endif
|
|
end subroutine potential_nh3_low
|
|
|
|
|
|
subroutine copy_2_lower_triangle(mat)
|
|
real(dp), intent(inout) :: mat(:, :)
|
|
integer :: m, n
|
|
! write lower triangle of matrix symmetrical
|
|
do n = 2, size(mat, 1)
|
|
do m = 1, n - 1
|
|
mat(n, m) = mat(m, n)
|
|
end do
|
|
end do
|
|
end subroutine copy_2_lower_triangle
|
|
|
|
|
|
end module
|