jnshuti
/
ANN-my-version
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

remove the model

This commit is contained in:
Jean Paul Nshuti 2026-06-18 09:21:49 +02:00
parent a7fef37ba5
commit 2c9c856bb4
16 changed files with 0 additions and 5399 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

149
src/model/Phase_correction.f90
View File

@ -1,149 +0,0 @@
module Phase_corr
! module for correcting the phase issue in the ab inition dipole data
use iso_fortran_env, only: dp => real64, idp => int32
use nncommons, only: ndiab
implicit none
!integer(idp), parameter:: ndiab = 4
contains
SUBROUTINE build_cluster(m,p)
! subroutine which check the sign of two succesive point
!and correct it if neccessary
! take dipole mx at point i and at point i+1
! update dipole at i+1
implicit none
Real(dp), intent(in):: m(ndiab,ndiab) ! at point i
Real(dp), intent(inout):: p(ndiab,ndiab) ! mx at i+1
!Real(dp), intent(inout):: q(qn)
!Real(dp):: prod2, prod3
! check first the sign of mx0(2,3) if it the same as sign of y(coord)
! q3 is ys
! q5 is yb
! check if the product of the transition element
! of two successive point there is no flip of sign
! check if state 1 need to be flipped
!if (needs_flip(M,P,1)) then
! write(6,*) "state 1 is flipped"
! call flip_e_state(P,1)
!endif
! check if e-state (2) need to be flipped
if (needs_flip(M,P,2)) then
!write(6,*) "state 2 is flipped"
call flip_e_state(P,2)
endif
! check if state (3) need to be flipped
!if (needs_flip(M,P,3)) then
! write(6,*) 'State 3 is flipped '
! call flip_e_state(P,3)
!endif
! check for the last state
!if (needs_flip(M,P,4)) then
! write(6,*) "state 4 is flipped"
! call flip_e_state(P,4)
!endif
END SUBROUTINE build_cluster
SUBROUTINE flip_e_state(A,O)
! the routine which flip the sign of the e-state element
!USE dim_parameter, only: ndiab
!USE accuracy_constants, only: idp, dp
implicit none
real(dp), intent(inout):: A(ndiab,ndiab) ! dipole matrix
integer(idp), intent(in):: O ! the index of the e-state to be flipped
real(dp):: P(ndiab,ndiab)
integer(idp):: i
!logical,save :: first_call =.true.
! the identity matrix will be used to create the transformation matrix P
P = 0.0d0
do i =1,ndiab
P(i,i)=1.0d0
enddo
! flipped the sign of the desired e state
P(o,o) = -1.0d0
! apply P on A
! A' = P**TAP
A=matmul(transpose(P),matmul(A,P))
! write the index of electronic state flipped
END SUBROUTINE
logical function needs_flip(m, p, s)
! Decide whether electronic state s should be flipped
implicit none
real(dp), intent(in) :: m(ndiab, ndiab), p(ndiab, ndiab)
integer(idp), intent(in) :: s
real(dp), parameter:: tol = 1.0d-9,eps_dip =1.0d-6
real(dp) :: m1,m2,p1,p2,m3,p3
needs_flip = .false.
select case (s)
case(1)
m1 = m(1,2) ; p1 = p(1,2)
m2 = m(1,3) ; p2 = p(1,3)
m3 = m(1,4) ; p3 = p(1,4)
if ((abs(m1) < eps_dip) .and. (abs(p1) < eps_dip)) then
m1 = m3
p1 = p3
else if ((abs(m2) < eps_dip) .and. (abs(p2) < eps_dip)) then
m2 = m3
p2 = p3
endif
case (2)
m1 = m(2,3) ; p1 = p(2,3)
m2 = m(2,4) ; p2 = P(2,4)
case (3)
m1 = m(2,3) ; p1 = p(2,3)
m2 = m(3,4) ; p2 = p(3,4)
case(4)
m1 = m(2,4) ; p1 = p(2,4)
m2 = m(3,4) ; p2 = p(3,4)
m3 = m(1,4) ; p3 = p(1,4)
if ((abs(m1) < eps_dip) .and. (abs(p1) < eps_dip)) then
m1 = m3
p1 = p3
else if ((abs(m2) < eps_dip) .and. (abs(p2) < eps_dip)) then
m2 = m3
p2 = p3
endif
case default
return
end select
if (abs(m1*p1) < tol .and. abs(m2*p2) < tol ) return
!if (abs(m2*p2) < tol ) return
!if (m1*p1 < 0.0_dp .and. m2*p2 < 0.0_dp) then
if (m1*p1 < 0.0_dp) then
needs_flip = .true.
end if
end function needs_flip
end module

1843
src/model/all.f90
View File

File diff suppressed because it is too large Load Diff

58
src/model/ctrans.f
View File

@ -1,58 +0,0 @@
! Author: jnshuti
! Created: 2025-10-24 16:44:03
! Last modified: 2026-01-12 11:20:55 jnshuti
subroutine cart2mode(qq)
use invariants_mod, only: invariants
use nn_params
use nncommons
implicit none
double precision, intent(inout):: qq(maxnin)
double precision :: q(12),inv(4)
double precision :: qm(6)
double precision, parameter :: m_N=25526.04237518618d0
double precision, parameter :: m_H=1837.152647439619d0
double precision, parameter :: Ang2Bohr=1.889726124565d0
double precision diff(12), ref_geom(12)
double precision q_mod(12)
! transformation matrix
double precision mode(12,6)
! final internal coord vector
include 'newmodes.f'
ref_geom(1:12)=0.d0
ref_geom( 4 )=1.022871078d0
ref_geom( 7 )=-.5d0*ref_geom(4)
!ref_geom( 8 )=sqrt(1.5d0)*ref_geom(4)
ref_geom( 8)=0.5d0*sqrt(3.0d0)*ref_geom(4)
ref_geom(10 )=-.5d0*ref_geom(4)
!ref_geom(11 )=-sqrt(1.5d0)*ref_geom(4)
ref_geom(11)=-0.5d0*sqrt(3.0d0)*ref_geom(4)
!massN=25526.04237518618
!massH=1837.152647439619
!
q(1:12)=qq(1:12)
q_mod = (1.0d0/Ang2Bohr)*q
! difference vector
diff=q_mod-ref_geom
! mass weighting
diff(1:3)=diff(1:3)*sqrt(m_N)
diff(4:12)=diff(4:12)*sqrt(m_H)
qm=matmul(transpose(mode),diff)
do i =1,6
if ( abs(qm(i)) .lt. 1.0d-6) qm(i) = 0.0d0
enddo
qm(3) = -qm(3)
call invariants(qm(1),qm(2),qm(3),qm(4),qm(5),inv)
qq(1:len_in)=inv(1:len_in)
qq(len_in+1:len_in+5)=qm(2:6)
end subroutine cart2mode

174
src/model/data_transform.f90
View File

@ -1,174 +0,0 @@
! Author: jnshuti
! Created: 2025-10-28 10:05:45
! Last modified: 2026-02-04 10:18:40 jnshuti
! Subroutine which transform the data in any desired format
module data_transf_mod
use iso_fortran_env, only: dp => real64, idp => int32
use nncommons, only: sets, ndata, inp_out,ndiab
use nn_params, only: maxpout,maxpats
use Phase_corr
implicit none
!integer(idp),parameter:: ndiab=4
private
public :: data_transform
contains
SUBROUTINE data_transform(pat_out,npat)
! IN: pat_out: the output patterns
! x1 the input pattern
! the aim it to transform pat_out in any disired form
implicit none
real(dp),intent(inout):: pat_out(maxpout,maxpats)
!real(dp),intent(inout):: ref_out(maxpout,maxpats)
integer(idp), intent(in):: npat ! the total number of point
! internal varibles
integer(idp) i,j,total ,ii,m,n
integer(idp), dimension(sets):: flipped_scan2, flipped_scan3
integer(idp), dimension(sets):: flipped_scan4
Real(dp),allocatable:: mux(:,:,:)
integer(idp):: scn_strt, scn_end,iref
logical, parameter:: flip =.true.
!write(6,*)"The total numbe of points used in data transform =",npat
! allocatable varables
allocate(mux(ndiab,ndiab,npat))
! define the scann need to be flipped
! list of all the scan need to flip the sign depending on the result of the molpro calculation
flipped_scan2 = 0; flipped_scan3 = 0; flipped_scan4 =0
flipped_scan3 = [1,3,4,5,6,7,8,9,10,11, &
12,13,14,15,16,17,18,19,20, &
21,22,24,25,54,55,56, &
58,59,60,61,62,63]
flipped_scan2 = [27,28,29,30,31,32,33,34,35, &
36,37,38,39,40,41,42,43,44, &
45,46,47,48,49,50,51,52,53, &
65,67,68,69,70,71, &
72,73,74]
flipped_scan4 = [23,64]
!flipped_scan2 = [77,80,118,120]
!flipped_scan3 = [56,111]
total=0
do i =1,sets
do j=1,ndata(i)
total = total+1
!pat_out(1:inp_out,total)=pat_out(1:inp_out,total)!-ref_en
! transform the y to matrix form
call Y2mat(pat_out(5:inp_out,total),Mux(:,:,total))
if (flip) then
if (j .gt. 2) call build_cluster(Mux(:,:,total-1),Mux(:,:,total))
endif
enddo
enddo
! flipe the required scan
if (flip) then
scn_end = 0
iref=0
do i =1,sets
if (ndata(i) .eq. 0) cycle
scn_strt = scn_end +1
scn_end = scn_end + ndata(i)
if (any(flipped_scan2 .eq. remap(i) )) then
do j = scn_strt,scn_end
call flip_e_state(Mux(:,:,j),2)
enddo
endif
if (any(flipped_scan3 .eq. remap(i) )) then
do j = scn_strt,scn_end
call flip_e_state(Mux(:,:,j),3)
enddo
endif
if (any(flipped_scan4 .eq. remap(i))) then
do j =scn_strt,scn_end
call flip_e_state(Mux(:,:,j),4)
enddo
endif
ii = 4
do m=1,ndiab
do n = m,ndiab
ii = ii +1
pat_out(ii,scn_strt:scn_end) = &
Mux(m,n,scn_strt:scn_end)
enddo
enddo
enddo
endif
deallocate(Mux)
END SUBROUTINE data_transform
subroutine Y2mat(Y,Mx)
IMPLICIT NONE
Real(dp), intent(in):: y(:)
Real(dp),intent(out):: Mx(ndiab,ndiab)
! internal variables
integer(idp):: ii,i,j
!ntot = ndiab*(ndiab+1)/2
!if (inp_out/2 .ne. ntot) then
! write(6,*)"output not equal ntot", inp_out,ntot
! stop
!endif
ii=1
do i=1,ndiab
do j=i,ndiab
! !mx
mx(i,j)=y(ii)
!
ii=ii+1
enddo
enddo
call copy_2_low(mx)
!call copy_2_low(my)
end subroutine
! ! subroutine which copy the upper matrix to lower
subroutine copy_2_low(m)
implicit none
real(dp), intent(inout) :: m(:,:)
integer(idp) :: i,j
! copy the lower part of the matrix to the upper part
do i=1,size(m,1)
do j=i,size(m,1)
m(j,i) = m(i,j)
enddo
enddo
end subroutine copy_2_low
integer function remap(i)
integer(idp),intent(in):: i
remap = i + 0
end function remap
end module data_transf_mod

411
src/model/dipole_model.f90
View File

@ -1,411 +0,0 @@
module diabmodel
use dip_param
use iso_fortran_env, only:dp => real64,idp => int32
use nn_params
use nncommons, only: len_in
implicit none
!integer(idp),parameter:: ndiab=4
!logical :: debug=.false.
private
public:: diab_x
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine diab_x(e,q,nn_out)
real(dp),intent(in) :: q(maxnin)
real(dp),intent(inout) :: nn_out(maxnout)
real(dp),intent(out) :: e(:,:)
integer(idp) id,i,j, ii
real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
real(dp), dimension(16) :: shift,scal
! real(dp),dimension(17):: shi,scalee
real(dp), parameter:: tol = 1.0d-9
!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)
id = 11
call init_dip_planar_data()
shift =1.0d0
scal = 1.0d0
ii = 5
do i = 28 ,np ! 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,*)" Dipole non-zero", ii-1
ss=xs**2+ys**2 ! totaly symmetric term
sb=xb**2+yb**2
v3_vec( 1) = xs*(xs**2-3*ys**2)
v3_vec( 2) = xb*(xb**2-3*yb**2)
v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys
v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb
v3_vec( 5) = ys*(3*xs**2-ys**2)
v3_vec( 6) = yb*(3*xb**2-yb**2)
v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys
v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb
e=0.0d0
!id=key !1
! V-term
! order 1
e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
id=id+1 !2
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
id=id+1 !3
e(4,4)=e(4,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
! order 2
id=id+1 !4
e(1,1)=e(1,1)+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)
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)
id=id+1 !6
e(4,4)=e(4,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)
! order 3
id=id+1 !7
e(1,1)=e(1,1)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb
id=id+1 !8
e(2,2)=e(2,2)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb
e(3,3)=e(3,3)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb
id=id+1 !9
e(4,4)=e(4,4)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb
! JAHN TELLER COUPLING W AND Z
! order 0
id=id+1 !10
e(2,2)=e(2,2)+p(pst(1,id))
e(3,3)=e(3,3)-p(pst(1,id))
! order 1
id=id+1 !11
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 !12
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)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb
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)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb)
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)
! order 3
id=id+1 !13
do i=1,4
j=i-1
e(2,2)=e(2,2)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i)
e(3,3)=e(3,3)-(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i)
e(2,3)=e(2,3)+(-p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i+4)
enddo
e(2,2)=e(2,2)+p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb
e(3,3)=e(3,3)-(p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb)
e(2,3)=e(2,3)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb
! PSEUDO JAHN TELLER
! A2 ground state coupled with E
! ###################################################
! ###################################################
! order 0
id=id+1 !14
e(1,2)=e(1,2)+b*p(pst(1,id))
! order 1
id=id+1 !15
e(1,2)=e(1,2)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
e(1,3)=e(1,3)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb)
! order 2
id=id+1 !16
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) + p(pst(1,id)+3)*(xs**2+ys**2))
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))
! order 3
id =id+1 ! 17
do i=1,4
e(1,2)=e(1,2)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
e(1,3)=e(1,3)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
enddo
!! THE COUPLING OF A2 WITH A1
!####################################################
!####################################################
! order 1
id=id+1 !18
e(1,4)=e(1,4)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
id=id+1 !19
e(1,4)=e(1,4)+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))
!!! THE COUPLING OF A1 WITH E
!!####################################################
!####################################################
! order 0
id=id+1 !20
e(2,4)=e(2,4)+p(pst(1,id))
! order 1
id=id+1 !21
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
! order 2
id=id+1 !22
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) +p(pst(1,id)+3)*(xs**2+ys**2)
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)
! order 3
id=id+1 !23
do i=1,4
e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
e(3,4)=e(3,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
enddo
!! End of the model
e(2,1)=e(1,2)
e(3,1)=e(1,3)
e(3,2)=e(2,3)
e(4,1)=e(1,4)
e(4,2)=e(2,4)
e(4,3)=e(3,4)
end subroutine diab_x
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! THE Y COMPONENT OF DIPOLE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!subroutine diab_y(e,q,nn_out)
! real(dp),intent(in) :: q(maxnin)
! real(dp),intent(inout) :: nn_out(maxnout)
! real(dp),intent(out) :: e(:,:)
! integer(idp) id,i,j, ii
! real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
! real(dp), dimension(23) :: shift,scal
! real(dp), parameter:: tol = 1.0d-9
!
! xs=q(5)
! ys=q(6)
! xb=q(7)
! yb=q(8)
! a=q(4)
! b=q(9)
!
! id = 11
!
! call init_dip_planar_data()
! scal = 1.0d-3
! shift =1.0d0
! ii = 1
! do i = 11 ,np ! 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
!
! id = 11
!
! ss=xs**2+ys**2 ! totaly symmetric term
! sb=xb**2+yb**2
!
! v3_vec( 1) = xs*(xs**2-3*ys**2)
! v3_vec( 2) = xb*(xb**2-3*yb**2)
! v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys
! v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb
! v3_vec( 5) = ys*(3*xs**2-ys**2)
! v3_vec( 6) = yb*(3*xb**2-yb**2)
! v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys
! v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb
!
!
!
!
!
! e=0.0d0
! ! V-term
! !id=key !1
! e(1,1)=e(1,1)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb
!
! id=id+1 !2
! e(2,2)=e(2,2)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb
! e(3,3)=e(3,3)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb
! id=id+1 !3
! e(4,4)=e(4,4)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb
! ! order 2
! id=id+1 !4
! e(1,1)=e(1,1)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) &
! -p(pst(1,id)+2)*(xs*yb+xb*ys)
! id=id+1 !5
! e(2,2)=e(2,2)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) &
! -p(pst(1,id)+2)*(xs*yb+xb*ys)
!
!
! e(3,3)=e(3,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)
! id=id+1 !6
! e(4,4)=e(4,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)
! ! order 3
! id=id+1 !7
! e(1,1)=e(1,1)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb
! id=id+1 !8
! e(2,2)=e(2,2)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb
! e(3,3)=e(3,3)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb
! id=id+1 !9
! e(4,4)=e(4,4)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb
!
! ! V- term + totally symmetric coord a
!
! ! JAHN TELLER COUPLING TERM
! ! order 0
! id=id+1 !10
! e(2,3)=e(2,3)+p(pst(1,id))
! ! order 1
!
! id=id+1 !11
! e(2,2)=e(2,2)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb
! e(3,3)=e(3,3)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb
! e(2,3)=e(2,3)-p(pst(1,id))*xs-p(pst(1,id)+1)*xb
! !id=id+1 !12
! ! order 2
! id=id+1 !12
! e(2,2)=e(2,2)+p(pst(1,id))*2*xs*ys+p(pst(1,id)+1)*2*xb*yb+p(pst(1,id)+2)*(xs*yb+xb*ys)
! e(3,3)=e(3,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)
! e(2,3)=e(2,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)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb
! ! order 3
! id=id+1 !13
! do i=1,4
! j=i-1
! e(2,2)=e(2,2)+(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4)
! e(3,3)=e(3,3)-(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4)
! e(2,3)=e(2,3)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i)
! enddo
! e(2,2)=e(2,2)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb
! e(3,3)=e(3,3)+p(pst(1,id)+8)*ys*ss+p(pst(1,id)+9)*yb*sb
! e(2,3)=e(2,3)-p(pst(1,id)+8)*xs*ss-p(pst(1,id)+1)*xb*sb
! ! PSEUDO JAHN TELLER
! ! ORDER 0
! ! THE COUPLING OF A2 GROUND STATE WITH E
! ! ###################################################
! ! ###################################################
! ! order 0
! id=id+1 !14
! e(1,3)=e(1,3)-b*(p(pst(1,id)))
! ! order 1
! id=id+1 !15
! e(1,2)=e(1,2)-b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb)
! e(1,3)=e(1,3)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
!
! ! order 2
! id=id+1 !16
! e(1,2)=e(1,2)+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))
! e(1,3)=e(1,3)+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) - p(pst(1,id)+3)*(xs**2+ys**2))
!
! ! order 3
! id = id+1 ! 17
! do i=1,4
! e(1,2)=e(1,2)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
! e(1,3)=e(1,3)-b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
! enddo
!
!
! ! THE COUPLING OF A2 WITH A1
! !####################################################
! !####################################################
! ! order 1
! id=id+1 !17
! e(1,4)=e(1,4)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb)
! ! order 2
! id=id+1 !18
! e(1,4)=e(1,4)-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 WITH E
! !####################################################
! !####################################################
! ! order 0
! id=id+1 !19
! e(3,4)=e(3,4)-p(pst(1,id))
! ! order 1
! id=id+1 !20
! e(2,4)=e(2,4)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb
! e(3,4)=e(3,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
! ! order 2
! id=id+1 !21
! e(2,4)=e(2,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)
! e(3,4)=e(3,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) - p(pst(1,id)+3)*(xs**2+ys**2)
!
! id =id+1 ! 23
! ! order 3
! do i=1,4
! e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
! e(3,4)=e(3,4)-(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
! enddo
!
! ! end of the model
! e(2,1)=e(1,2)
! e(3,1)=e(1,3)
! e(3,2)=e(2,3)
! e(4,1)=e(1,4)
! e(4,2)=e(2,4)
! e(4,3)=e(3,4)
! ! change Ex and Ey
!end subroutine diab_y
! potential nh3
! parth
!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 diabmodel

914
src/model/fit_genetic.f90
View File

@ -1,914 +0,0 @@
Module dip_param
IMPLICIT NONE
Integer,parameter :: np=101
Double precision :: p(101)
integer :: pst(2,400)
! Non zero parameters are 17
!$omp threadprivate(p,pst)
contains
SUBROUTINE init_dip_planar_data()
implicit none
p( 1)= -0.561433980D+02
p( 2)= 0.000000000D+00
p( 3)= 0.000000000D+00
p( 4)= -0.558830090D+02
p( 5)= 0.000000000D+00
p( 6)= 0.000000000D+00
p( 7)= -0.557220030D+02
p( 8)= 0.000000000D+00
p( 9)= 0.000000000D+00
p( 10)= 0.000000000D+00
p( 11)= 0.334452655D-02
p( 12)= 0.000000000D+00
p( 13)= 0.000000000D+00
p( 14)= 0.000000000D+00
p( 15)= 0.000000000D+00
p( 16)= 0.000000000D+00
p( 17)= 0.000000000D+00
p( 18)= 0.000000000D+00
p( 19)= 0.000000000D+00
p( 20)= 0.000000000D+00
p( 21)= -0.532590808D-03
p( 22)= 0.000000000D+00
p( 23)= 0.000000000D+00
p( 24)= 0.000000000D+00
p( 25)= 0.000000000D+00
p( 26)= 0.000000000D+00
p( 27)= 0.000000000D+00
p( 28)= -0.226116045D-01
p( 29)= 0.000000000D+00
p( 30)= -0.363599361D-01
p( 31)= 0.000000000D+00
p( 32)= 0.962100323D-01
p( 33)= 0.000000000D+00
p( 34)= 0.111336127D-03
p( 35)= 0.000000000D+00
p( 36)= 0.000000000D+00
p( 37)= -0.433685770D-04
p( 38)= 0.000000000D+00
p( 39)= 0.000000000D+00
p( 40)= -0.548870572D-02
p( 41)= 0.000000000D+00
p( 42)= 0.000000000D+00
p( 43)= 0.000000000D+00
p( 44)= 0.000000000D+00
p( 45)= 0.000000000D+00
p( 46)= 0.000000000D+00
p( 47)= 0.000000000D+00
p( 48)= 0.000000000D+00
p( 49)= 0.217598103D+00
p( 50)= -0.113909871D-01
p( 51)= 0.000000000D+00
p( 52)= -0.123184224D-03
p( 53)= 0.000000000D+00
p( 54)= 0.000000000D+00
p( 55)= 0.000000000D+00
p( 56)= 0.000000000D+00
p( 57)= 0.000000000D+00
p( 58)= 0.000000000D+00
p( 59)= 0.000000000D+00
p( 60)= 0.000000000D+00
p( 61)= 0.000000000D+00
p( 62)= 0.000000000D+00
p( 63)= 0.000000000D+00
p( 64)= 0.000000000D+00
p( 65)= 0.000000000D+00
p( 66)= 0.000000000D+00
p( 67)= 0.000000000D+00
p( 68)= 0.000000000D+00
p( 69)= 0.000000000D+00
p( 70)= 0.000000000D+00
p( 71)= 0.000000000D+00
p( 72)= 0.000000000D+00
p( 73)= 0.000000000D+00
p( 74)= 0.000000000D+00
p( 75)= 0.000000000D+00
p( 76)= 0.000000000D+00
p( 77)= 0.000000000D+00
p( 78)= 0.000000000D+00
p( 79)= 0.000000000D+00
p( 80)= 0.000000000D+00
p( 81)= 0.000000000D+00
p( 82)= 0.000000000D+00
p( 83)= 0.000000000D+00
p( 84)= 0.000000000D+00
p( 85)= 0.000000000D+00
p( 86)= 0.000000000D+00
p( 87)= -0.567206500D-02
p( 88)= -0.613490432D-02
p( 89)= 0.000000000D+00
p( 90)= 0.303022045D-03
p( 91)= 0.000000000D+00
p( 92)= 0.000000000D+00
p( 93)= 0.000000000D+00
p( 94)= 0.000000000D+00
p( 95)= 0.000000000D+00
p( 96)= 0.000000000D+00
p( 97)= 0.000000000D+00
p( 98)= 0.000000000D+00
p( 99)= 0.000000000D+00
p(100)= 0.000000000D+00
p(101)= 0.000000000D+00
pst(1, 1)= 1
pst(2, 1)= 3
pst(1, 2)= 4
pst(2, 2)= 3
pst(1, 3)= 7
pst(2, 3)= 4
pst(1, 4)= 11
pst(2, 4)= 2
pst(1, 5)= 13
pst(2, 5)= 3
pst(1, 6)= 16
pst(2, 6)= 2
pst(1, 7)= 18
pst(2, 7)= 3
pst(1, 8)= 21
pst(2, 8)= 3
pst(1, 9)= 24
pst(2, 9)= 3
pst(1, 10)= 27
pst(2, 10)= 1
pst(1, 11)= 28
pst(2, 11)= 2
pst(1, 12)= 30
pst(2, 12)= 2
pst(1, 13)= 32
pst(2, 13)= 2
pst(1, 14)= 34
pst(2, 14)= 3
pst(1, 15)= 37
pst(2, 15)= 3
pst(1, 16)= 40
pst(2, 16)= 3
pst(1, 17)= 43
pst(2, 17)= 2
pst(1, 18)= 45
pst(2, 18)= 2
pst(1, 19)= 47
pst(2, 19)= 2
pst(1, 20)= 49
pst(2, 20)= 1
pst(1, 21)= 50
pst(2, 21)= 2
pst(1, 22)= 52
pst(2, 22)= 5
pst(1, 23)= 57
pst(2, 23)= 10
pst(1, 24)= 67
pst(2, 24)= 1
pst(1, 25)= 68
pst(2, 25)= 2
pst(1, 26)= 70
pst(2, 26)= 4
pst(1, 27)= 74
pst(2, 27)= 8
pst(1, 28)= 82
pst(2, 28)= 2
pst(1, 29)= 84
pst(2, 29)= 3
pst(1, 30)= 87
pst(2, 30)= 1
pst(1, 31)= 88
pst(2, 31)= 2
pst(1, 32)= 90
pst(2, 32)= 4
pst(1, 33)= 94
pst(2, 33)= 8
pst(1, 34)= 0
pst(2, 34)= 0
pst(1, 35)= 0
pst(2, 35)= 0
pst(1, 36)= 0
pst(2, 36)= 0
pst(1, 37)= 0
pst(2, 37)= 0
pst(1, 38)= 0
pst(2, 38)= 0
pst(1, 39)= 0
pst(2, 39)= 0
pst(1, 40)= 0
pst(2, 40)= 0
pst(1, 41)= 0
pst(2, 41)= 0
pst(1, 42)= 0
pst(2, 42)= 0
pst(1, 43)= 0
pst(2, 43)= 0
pst(1, 44)= 0
pst(2, 44)= 0
pst(1, 45)= 0
pst(2, 45)= 0
pst(1, 46)= 0
pst(2, 46)= 0
pst(1, 47)= 0
pst(2, 47)= 0
pst(1, 48)= 0
pst(2, 48)= 0
pst(1, 49)= 0
pst(2, 49)= 0
pst(1, 50)= 0
pst(2, 50)= 0
pst(1, 51)= 0
pst(2, 51)= 0
pst(1, 52)= 0
pst(2, 52)= 0
pst(1, 53)= 0
pst(2, 53)= 0
pst(1, 54)= 0
pst(2, 54)= 0
pst(1, 55)= 0
pst(2, 55)= 0
pst(1, 56)= 0
pst(2, 56)= 0
pst(1, 57)= 0
pst(2, 57)= 0
pst(1, 58)= 0
pst(2, 58)= 0
pst(1, 59)= 0
pst(2, 59)= 0
pst(1, 60)= 0
pst(2, 60)= 0
pst(1, 61)= 0
pst(2, 61)= 0
pst(1, 62)= 0
pst(2, 62)= 0
pst(1, 63)= 0
pst(2, 63)= 0
pst(1, 64)= 0
pst(2, 64)= 0
pst(1, 65)= 0
pst(2, 65)= 0
pst(1, 66)= 0
pst(2, 66)= 0
pst(1, 67)= 0
pst(2, 67)= 0
pst(1, 68)= 0
pst(2, 68)= 0
pst(1, 69)= 0
pst(2, 69)= 0
pst(1, 70)= 0
pst(2, 70)= 0
pst(1, 71)= 0
pst(2, 71)= 0
pst(1, 72)= 0
pst(2, 72)= 0
pst(1, 73)= 0
pst(2, 73)= 0
pst(1, 74)= 0
pst(2, 74)= 0
pst(1, 75)= 0
pst(2, 75)= 0
pst(1, 76)= 0
pst(2, 76)= 0
pst(1, 77)= 0
pst(2, 77)= 0
pst(1, 78)= 0
pst(2, 78)= 0
pst(1, 79)= 0
pst(2, 79)= 0
pst(1, 80)= 0
pst(2, 80)= 0
pst(1, 81)= 0
pst(2, 81)= 0
pst(1, 82)= 0
pst(2, 82)= 0
pst(1, 83)= 0
pst(2, 83)= 0
pst(1, 84)= 0
pst(2, 84)= 0
pst(1, 85)= 0
pst(2, 85)= 0
pst(1, 86)= 0
pst(2, 86)= 0
pst(1, 87)= 0
pst(2, 87)= 0
pst(1, 88)= 0
pst(2, 88)= 0
pst(1, 89)= 0
pst(2, 89)= 0
pst(1, 90)= 0
pst(2, 90)= 0
pst(1, 91)= 0
pst(2, 91)= 0
pst(1, 92)= 0
pst(2, 92)= 0
pst(1, 93)= 0
pst(2, 93)= 0
pst(1, 94)= 0
pst(2, 94)= 0
pst(1, 95)= 0
pst(2, 95)= 0
pst(1, 96)= 0
pst(2, 96)= 0
pst(1, 97)= 0
pst(2, 97)= 0
pst(1, 98)= 0
pst(2, 98)= 0
pst(1, 99)= 0
pst(2, 99)= 0
pst(1,100)= 0
pst(2,100)= 0
pst(1,101)= 0
pst(2,101)= 0
pst(1,102)= 0
pst(2,102)= 0
pst(1,103)= 0
pst(2,103)= 0
pst(1,104)= 0
pst(2,104)= 0
pst(1,105)= 0
pst(2,105)= 0
pst(1,106)= 0
pst(2,106)= 0
pst(1,107)= 0
pst(2,107)= 0
pst(1,108)= 0
pst(2,108)= 0
pst(1,109)= 0
pst(2,109)= 0
pst(1,110)= 0
pst(2,110)= 0
pst(1,111)= 0
pst(2,111)= 0
pst(1,112)= 0
pst(2,112)= 0
pst(1,113)= 0
pst(2,113)= 0
pst(1,114)= 0
pst(2,114)= 0
pst(1,115)= 0
pst(2,115)= 0
pst(1,116)= 0
pst(2,116)= 0
pst(1,117)= 0
pst(2,117)= 0
pst(1,118)= 0
pst(2,118)= 0
pst(1,119)= 0
pst(2,119)= 0
pst(1,120)= 0
pst(2,120)= 0
pst(1,121)= 0
pst(2,121)= 0
pst(1,122)= 0
pst(2,122)= 0
pst(1,123)= 0
pst(2,123)= 0
pst(1,124)= 0
pst(2,124)= 0
pst(1,125)= 0
pst(2,125)= 0
pst(1,126)= 0
pst(2,126)= 0
pst(1,127)= 0
pst(2,127)= 0
pst(1,128)= 0
pst(2,128)= 0
pst(1,129)= 0
pst(2,129)= 0
pst(1,130)= 0
pst(2,130)= 0
pst(1,131)= 0
pst(2,131)= 0
pst(1,132)= 0
pst(2,132)= 0
pst(1,133)= 0
pst(2,133)= 0
pst(1,134)= 0
pst(2,134)= 0
pst(1,135)= 0
pst(2,135)= 0
pst(1,136)= 0
pst(2,136)= 0
pst(1,137)= 0
pst(2,137)= 0
pst(1,138)= 0
pst(2,138)= 0
pst(1,139)= 0
pst(2,139)= 0
pst(1,140)= 0
pst(2,140)= 0
pst(1,141)= 0
pst(2,141)= 0
pst(1,142)= 0
pst(2,142)= 0
pst(1,143)= 0
pst(2,143)= 0
pst(1,144)= 0
pst(2,144)= 0
pst(1,145)= 0
pst(2,145)= 0
pst(1,146)= 0
pst(2,146)= 0
pst(1,147)= 0
pst(2,147)= 0
pst(1,148)= 0
pst(2,148)= 0
pst(1,149)= 0
pst(2,149)= 0
pst(1,150)= 0
pst(2,150)= 0
pst(1,151)= 0
pst(2,151)= 0
pst(1,152)= 0
pst(2,152)= 0
pst(1,153)= 0
pst(2,153)= 0
pst(1,154)= 0
pst(2,154)= 0
pst(1,155)= 0
pst(2,155)= 0
pst(1,156)= 0
pst(2,156)= 0
pst(1,157)= 0
pst(2,157)= 0
pst(1,158)= 0
pst(2,158)= 0
pst(1,159)= 0
pst(2,159)= 0
pst(1,160)= 0
pst(2,160)= 0
pst(1,161)= 0
pst(2,161)= 0
pst(1,162)= 0
pst(2,162)= 0
pst(1,163)= 0
pst(2,163)= 0
pst(1,164)= 0
pst(2,164)= 0
pst(1,165)= 0
pst(2,165)= 0
pst(1,166)= 0
pst(2,166)= 0
pst(1,167)= 0
pst(2,167)= 0
pst(1,168)= 0
pst(2,168)= 0
pst(1,169)= 0
pst(2,169)= 0
pst(1,170)= 0
pst(2,170)= 0
pst(1,171)= 0
pst(2,171)= 0
pst(1,172)= 0
pst(2,172)= 0
pst(1,173)= 0
pst(2,173)= 0
pst(1,174)= 0
pst(2,174)= 0
pst(1,175)= 0
pst(2,175)= 0
pst(1,176)= 0
pst(2,176)= 0
pst(1,177)= 0
pst(2,177)= 0
pst(1,178)= 0
pst(2,178)= 0
pst(1,179)= 0
pst(2,179)= 0
pst(1,180)= 0
pst(2,180)= 0
pst(1,181)= 0
pst(2,181)= 0
pst(1,182)= 0
pst(2,182)= 0
pst(1,183)= 0
pst(2,183)= 0
pst(1,184)= 0
pst(2,184)= 0
pst(1,185)= 0
pst(2,185)= 0
pst(1,186)= 0
pst(2,186)= 0
pst(1,187)= 0
pst(2,187)= 0
pst(1,188)= 0
pst(2,188)= 0
pst(1,189)= 0
pst(2,189)= 0
pst(1,190)= 0
pst(2,190)= 0
pst(1,191)= 0
pst(2,191)= 0
pst(1,192)= 0
pst(2,192)= 0
pst(1,193)= 0
pst(2,193)= 0
pst(1,194)= 0
pst(2,194)= 0
pst(1,195)= 0
pst(2,195)= 0
pst(1,196)= 0
pst(2,196)= 0
pst(1,197)= 0
pst(2,197)= 0
pst(1,198)= 0
pst(2,198)= 0
pst(1,199)= 0
pst(2,199)= 0
pst(1,200)= 0
pst(2,200)= 0
pst(1,201)= 0
pst(2,201)= 0
pst(1,202)= 0
pst(2,202)= 0
pst(1,203)= 0
pst(2,203)= 0
pst(1,204)= 0
pst(2,204)= 0
pst(1,205)= 0
pst(2,205)= 0
pst(1,206)= 0
pst(2,206)= 0
pst(1,207)= 0
pst(2,207)= 0
pst(1,208)= 0
pst(2,208)= 0
pst(1,209)= 0
pst(2,209)= 0
pst(1,210)= 0
pst(2,210)= 0
pst(1,211)= 0
pst(2,211)= 0
pst(1,212)= 0
pst(2,212)= 0
pst(1,213)= 0
pst(2,213)= 0
pst(1,214)= 0
pst(2,214)= 0
pst(1,215)= 0
pst(2,215)= 0
pst(1,216)= 0
pst(2,216)= 0
pst(1,217)= 0
pst(2,217)= 0
pst(1,218)= 0
pst(2,218)= 0
pst(1,219)= 0
pst(2,219)= 0
pst(1,220)= 0
pst(2,220)= 0
pst(1,221)= 0
pst(2,221)= 0
pst(1,222)= 0
pst(2,222)= 0
pst(1,223)= 0
pst(2,223)= 0
pst(1,224)= 0
pst(2,224)= 0
pst(1,225)= 0
pst(2,225)= 0
pst(1,226)= 0
pst(2,226)= 0
pst(1,227)= 0
pst(2,227)= 0
pst(1,228)= 0
pst(2,228)= 0
pst(1,229)= 0
pst(2,229)= 0
pst(1,230)= 0
pst(2,230)= 0
pst(1,231)= 0
pst(2,231)= 0
pst(1,232)= 0
pst(2,232)= 0
pst(1,233)= 0
pst(2,233)= 0
pst(1,234)= 0
pst(2,234)= 0
pst(1,235)= 0
pst(2,235)= 0
pst(1,236)= 0
pst(2,236)= 0
pst(1,237)= 0
pst(2,237)= 0
pst(1,238)= 0
pst(2,238)= 0
pst(1,239)= 0
pst(2,239)= 0
pst(1,240)= 0
pst(2,240)= 0
pst(1,241)= 0
pst(2,241)= 0
pst(1,242)= 0
pst(2,242)= 0
pst(1,243)= 0
pst(2,243)= 0
pst(1,244)= 0
pst(2,244)= 0
pst(1,245)= 0
pst(2,245)= 0
pst(1,246)= 0
pst(2,246)= 0
pst(1,247)= 0
pst(2,247)= 0
pst(1,248)= 0
pst(2,248)= 0
pst(1,249)= 0
pst(2,249)= 0
pst(1,250)= 0
pst(2,250)= 0
pst(1,251)= 0
pst(2,251)= 0
pst(1,252)= 0
pst(2,252)= 0
pst(1,253)= 0
pst(2,253)= 0
pst(1,254)= 0
pst(2,254)= 0
pst(1,255)= 0
pst(2,255)= 0
pst(1,256)= 0
pst(2,256)= 0
pst(1,257)= 0
pst(2,257)= 0
pst(1,258)= 0
pst(2,258)= 0
pst(1,259)= 0
pst(2,259)= 0
pst(1,260)= 0
pst(2,260)= 0
pst(1,261)= 0
pst(2,261)= 0
pst(1,262)= 0
pst(2,262)= 0
pst(1,263)= 0
pst(2,263)= 0
pst(1,264)= 0
pst(2,264)= 0
pst(1,265)= 0
pst(2,265)= 0
pst(1,266)= 0
pst(2,266)= 0
pst(1,267)= 0
pst(2,267)= 0
pst(1,268)= 0
pst(2,268)= 0
pst(1,269)= 0
pst(2,269)= 0
pst(1,270)= 0
pst(2,270)= 0
pst(1,271)= 0
pst(2,271)= 0
pst(1,272)= 0
pst(2,272)= 0
pst(1,273)= 0
pst(2,273)= 0
pst(1,274)= 0
pst(2,274)= 0
pst(1,275)= 0
pst(2,275)= 0
pst(1,276)= 0
pst(2,276)= 0
pst(1,277)= 0
pst(2,277)= 0
pst(1,278)= 0
pst(2,278)= 0
pst(1,279)= 0
pst(2,279)= 0
pst(1,280)= 0
pst(2,280)= 0
pst(1,281)= 0
pst(2,281)= 0
pst(1,282)= 0
pst(2,282)= 0
pst(1,283)= 0
pst(2,283)= 0
pst(1,284)= 0
pst(2,284)= 0
pst(1,285)= 0
pst(2,285)= 0
pst(1,286)= 0
pst(2,286)= 0
pst(1,287)= 0
pst(2,287)= 0
pst(1,288)= 0
pst(2,288)= 0
pst(1,289)= 0
pst(2,289)= 0
pst(1,290)= 0
pst(2,290)= 0
pst(1,291)= 0
pst(2,291)= 0
pst(1,292)= 0
pst(2,292)= 0
pst(1,293)= 0
pst(2,293)= 0
pst(1,294)= 0
pst(2,294)= 0
pst(1,295)= 0
pst(2,295)= 0
pst(1,296)= 0
pst(2,296)= 0
pst(1,297)= 0
pst(2,297)= 0
pst(1,298)= 0
pst(2,298)= 0
pst(1,299)= 0
pst(2,299)= 0
pst(1,300)= 0
pst(2,300)= 0
pst(1,301)= 0
pst(2,301)= 0
pst(1,302)= 0
pst(2,302)= 0
pst(1,303)= 0
pst(2,303)= 0
pst(1,304)= 0
pst(2,304)= 0
pst(1,305)= 0
pst(2,305)= 0
pst(1,306)= 0
pst(2,306)= 0
pst(1,307)= 0
pst(2,307)= 0
pst(1,308)= 0
pst(2,308)= 0
pst(1,309)= 0
pst(2,309)= 0
pst(1,310)= 0
pst(2,310)= 0
pst(1,311)= 0
pst(2,311)= 0
pst(1,312)= 0
pst(2,312)= 0
pst(1,313)= 0
pst(2,313)= 0
pst(1,314)= 0
pst(2,314)= 0
pst(1,315)= 0
pst(2,315)= 0
pst(1,316)= 0
pst(2,316)= 0
pst(1,317)= 0
pst(2,317)= 0
pst(1,318)= 0
pst(2,318)= 0
pst(1,319)= 0
pst(2,319)= 0
pst(1,320)= 0
pst(2,320)= 0
pst(1,321)= 0
pst(2,321)= 0
pst(1,322)= 0
pst(2,322)= 0
pst(1,323)= 0
pst(2,323)= 0
pst(1,324)= 0
pst(2,324)= 0
pst(1,325)= 0
pst(2,325)= 0
pst(1,326)= 0
pst(2,326)= 0
pst(1,327)= 0
pst(2,327)= 0
pst(1,328)= 0
pst(2,328)= 0
pst(1,329)= 0
pst(2,329)= 0
pst(1,330)= 0
pst(2,330)= 0
pst(1,331)= 0
pst(2,331)= 0
pst(1,332)= 0
pst(2,332)= 0
pst(1,333)= 0
pst(2,333)= 0
pst(1,334)= 0
pst(2,334)= 0
pst(1,335)= 0
pst(2,335)= 0
pst(1,336)= 0
pst(2,336)= 0
pst(1,337)= 0
pst(2,337)= 0
pst(1,338)= 0
pst(2,338)= 0
pst(1,339)= 0
pst(2,339)= 0
pst(1,340)= 0
pst(2,340)= 0
pst(1,341)= 0
pst(2,341)= 0
pst(1,342)= 0
pst(2,342)= 0
pst(1,343)= 0
pst(2,343)= 0
pst(1,344)= 0
pst(2,344)= 0
pst(1,345)= 0
pst(2,345)= 0
pst(1,346)= 0
pst(2,346)= 0
pst(1,347)= 0
pst(2,347)= 0
pst(1,348)= 0
pst(2,348)= 0
pst(1,349)= 0
pst(2,349)= 0
pst(1,350)= 0
pst(2,350)= 0
pst(1,351)= 0
pst(2,351)= 0
pst(1,352)= 0
pst(2,352)= 0
pst(1,353)= 0
pst(2,353)= 0
pst(1,354)= 0
pst(2,354)= 0
pst(1,355)= 0
pst(2,355)= 0
pst(1,356)= 0
pst(2,356)= 0
pst(1,357)= 0
pst(2,357)= 0
pst(1,358)= 0
pst(2,358)= 0
pst(1,359)= 0
pst(2,359)= 0
pst(1,360)= 0
pst(2,360)= 0
pst(1,361)= 0
pst(2,361)= 0
pst(1,362)= 0
pst(2,362)= 0
pst(1,363)= 0
pst(2,363)= 0
pst(1,364)= 0
pst(2,364)= 0
pst(1,365)= 0
pst(2,365)= 0
pst(1,366)= 0
pst(2,366)= 0
pst(1,367)= 0
pst(2,367)= 0
pst(1,368)= 0
pst(2,368)= 0
pst(1,369)= 0
pst(2,369)= 0
pst(1,370)= 0
pst(2,370)= 0
pst(1,371)= 0
pst(2,371)= 0
pst(1,372)= 0
pst(2,372)= 0
pst(1,373)= 0
pst(2,373)= 0
pst(1,374)= 0
pst(2,374)= 0
pst(1,375)= 0
pst(2,375)= 0
pst(1,376)= 0
pst(2,376)= 0
pst(1,377)= 0
pst(2,377)= 0
pst(1,378)= 0
pst(2,378)= 0
pst(1,379)= 0
pst(2,379)= 0
pst(1,380)= 0
pst(2,380)= 0
pst(1,381)= 0
pst(2,381)= 0
pst(1,382)= 0
pst(2,382)= 0
pst(1,383)= 0
pst(2,383)= 0
pst(1,384)= 0
pst(2,384)= 0
pst(1,385)= 0
pst(2,385)= 0
pst(1,386)= 0
pst(2,386)= 0
pst(1,387)= 0
pst(2,387)= 0
pst(1,388)= 0
pst(2,388)= 0
pst(1,389)= 0
pst(2,389)= 0
pst(1,390)= 0
pst(2,390)= 0
pst(1,391)= 0
pst(2,391)= 0
pst(1,392)= 0
pst(2,392)= 0
pst(1,393)= 0
pst(2,393)= 0
pst(1,394)= 0
pst(2,394)= 0
pst(1,395)= 0
pst(2,395)= 0
pst(1,396)= 0
pst(2,396)= 0
pst(1,397)= 0
pst(2,397)= 0
pst(1,398)= 0
pst(2,398)= 0
pst(1,399)= 0
pst(2,399)= 0
pst(1,400)= 0
pst(2,400)= 0
End SUBROUTINE init_dip_planar_data
End Module dip_param

914
src/model/fit_param_dip_pes.f90
View File

@ -1,914 +0,0 @@
Module dip_param
IMPLICIT NONE
Integer,parameter :: np=101
Double precision :: p(101)
integer :: pst(2,400)
! Non zero parameters are 16
!$omp threadprivate(p,pst)
contains
SUBROUTINE init_dip_planar_data()
implicit none
p( 1)= -0.561433980D+02
p( 2)= 0.000000000D+00
p( 3)= 0.000000000D+00
p( 4)= -0.558830090D+02
p( 5)= 0.000000000D+00
p( 6)= 0.000000000D+00
p( 7)= -0.557220030D+02
p( 8)= 0.000000000D+00
p( 9)= 0.000000000D+00
p( 10)= 0.000000000D+00
p( 11)= 0.333624014D-02
p( 12)= 0.000000000D+00
p( 13)= 0.000000000D+00
p( 14)= 0.000000000D+00
p( 15)= 0.000000000D+00
p( 16)= 0.000000000D+00
p( 17)= 0.000000000D+00
p( 18)= 0.000000000D+00
p( 19)= 0.000000000D+00
p( 20)= 0.000000000D+00
p( 21)= -0.532590808D-03
p( 22)= 0.000000000D+00
p( 23)= 0.000000000D+00
p( 24)= 0.000000000D+00
p( 25)= 0.000000000D+00
p( 26)= 0.000000000D+00
p( 27)= 0.000000000D+00
p( 28)= -0.226073502D-01
p( 29)= 0.000000000D+00
p( 30)= -0.363395355D-01
p( 31)= 0.000000000D+00
p( 32)= 0.939930347D-01
p( 33)= 0.000000000D+00
p( 34)= 0.111364489D-03
p( 35)= 0.000000000D+00
p( 36)= 0.000000000D+00
p( 37)= 0.000000000D+00
p( 38)= 0.000000000D+00
p( 39)= 0.000000000D+00
p( 40)= -0.545384229D-02
p( 41)= 0.000000000D+00
p( 42)= 0.000000000D+00
p( 43)= 0.000000000D+00
p( 44)= 0.000000000D+00
p( 45)= 0.000000000D+00
p( 46)= 0.000000000D+00
p( 47)= 0.000000000D+00
p( 48)= 0.000000000D+00
p( 49)= 0.217598103D+00
p( 50)= -0.113772250D-01
p( 51)= 0.000000000D+00
p( 52)= -0.122160662D-03
p( 53)= 0.000000000D+00
p( 54)= 0.000000000D+00
p( 55)= 0.000000000D+00
p( 56)= 0.000000000D+00
p( 57)= 0.000000000D+00
p( 58)= 0.000000000D+00
p( 59)= 0.000000000D+00
p( 60)= 0.000000000D+00
p( 61)= 0.000000000D+00
p( 62)= 0.000000000D+00
p( 63)= 0.000000000D+00
p( 64)= 0.000000000D+00
p( 65)= 0.000000000D+00
p( 66)= 0.000000000D+00
p( 67)= 0.000000000D+00
p( 68)= 0.000000000D+00
p( 69)= 0.000000000D+00
p( 70)= 0.000000000D+00
p( 71)= 0.000000000D+00
p( 72)= 0.000000000D+00
p( 73)= 0.000000000D+00
p( 74)= 0.000000000D+00
p( 75)= 0.000000000D+00
p( 76)= 0.000000000D+00
p( 77)= 0.000000000D+00
p( 78)= 0.000000000D+00
p( 79)= 0.000000000D+00
p( 80)= 0.000000000D+00
p( 81)= 0.000000000D+00
p( 82)= 0.000000000D+00
p( 83)= 0.000000000D+00
p( 84)= 0.000000000D+00
p( 85)= 0.000000000D+00
p( 86)= 0.000000000D+00
p( 87)= -0.567206500D-02
p( 88)= -0.565873672D-02
p( 89)= 0.000000000D+00
p( 90)= 0.336524505D-03
p( 91)= 0.000000000D+00
p( 92)= 0.000000000D+00
p( 93)= 0.000000000D+00
p( 94)= 0.000000000D+00
p( 95)= 0.000000000D+00
p( 96)= 0.000000000D+00
p( 97)= 0.000000000D+00
p( 98)= 0.000000000D+00
p( 99)= 0.000000000D+00
p(100)= 0.000000000D+00
p(101)= 0.000000000D+00
pst(1, 1)= 1
pst(2, 1)= 3
pst(1, 2)= 4
pst(2, 2)= 3
pst(1, 3)= 7
pst(2, 3)= 4
pst(1, 4)= 11
pst(2, 4)= 2
pst(1, 5)= 13
pst(2, 5)= 3
pst(1, 6)= 16
pst(2, 6)= 2
pst(1, 7)= 18
pst(2, 7)= 3
pst(1, 8)= 21
pst(2, 8)= 3
pst(1, 9)= 24
pst(2, 9)= 3
pst(1, 10)= 27
pst(2, 10)= 1
pst(1, 11)= 28
pst(2, 11)= 2
pst(1, 12)= 30
pst(2, 12)= 2
pst(1, 13)= 32
pst(2, 13)= 2
pst(1, 14)= 34
pst(2, 14)= 3
pst(1, 15)= 37
pst(2, 15)= 3
pst(1, 16)= 40
pst(2, 16)= 3
pst(1, 17)= 43
pst(2, 17)= 2
pst(1, 18)= 45
pst(2, 18)= 2
pst(1, 19)= 47
pst(2, 19)= 2
pst(1, 20)= 49
pst(2, 20)= 1
pst(1, 21)= 50
pst(2, 21)= 2
pst(1, 22)= 52
pst(2, 22)= 5
pst(1, 23)= 57
pst(2, 23)= 10
pst(1, 24)= 67
pst(2, 24)= 1
pst(1, 25)= 68
pst(2, 25)= 2
pst(1, 26)= 70
pst(2, 26)= 4
pst(1, 27)= 74
pst(2, 27)= 8
pst(1, 28)= 82
pst(2, 28)= 2
pst(1, 29)= 84
pst(2, 29)= 3
pst(1, 30)= 87
pst(2, 30)= 1
pst(1, 31)= 88
pst(2, 31)= 2
pst(1, 32)= 90
pst(2, 32)= 4
pst(1, 33)= 94
pst(2, 33)= 8
pst(1, 34)= 0
pst(2, 34)= 0
pst(1, 35)= 0
pst(2, 35)= 0
pst(1, 36)= 0
pst(2, 36)= 0
pst(1, 37)= 0
pst(2, 37)= 0
pst(1, 38)= 0
pst(2, 38)= 0
pst(1, 39)= 0
pst(2, 39)= 0
pst(1, 40)= 0
pst(2, 40)= 0
pst(1, 41)= 0
pst(2, 41)= 0
pst(1, 42)= 0
pst(2, 42)= 0
pst(1, 43)= 0
pst(2, 43)= 0
pst(1, 44)= 0
pst(2, 44)= 0
pst(1, 45)= 0
pst(2, 45)= 0
pst(1, 46)= 0
pst(2, 46)= 0
pst(1, 47)= 0
pst(2, 47)= 0
pst(1, 48)= 0
pst(2, 48)= 0
pst(1, 49)= 0
pst(2, 49)= 0
pst(1, 50)= 0
pst(2, 50)= 0
pst(1, 51)= 0
pst(2, 51)= 0
pst(1, 52)= 0
pst(2, 52)= 0
pst(1, 53)= 0
pst(2, 53)= 0
pst(1, 54)= 0
pst(2, 54)= 0
pst(1, 55)= 0
pst(2, 55)= 0
pst(1, 56)= 0
pst(2, 56)= 0
pst(1, 57)= 0
pst(2, 57)= 0
pst(1, 58)= 0
pst(2, 58)= 0
pst(1, 59)= 0
pst(2, 59)= 0
pst(1, 60)= 0
pst(2, 60)= 0
pst(1, 61)= 0
pst(2, 61)= 0
pst(1, 62)= 0
pst(2, 62)= 0
pst(1, 63)= 0
pst(2, 63)= 0
pst(1, 64)= 0
pst(2, 64)= 0
pst(1, 65)= 0
pst(2, 65)= 0
pst(1, 66)= 0
pst(2, 66)= 0
pst(1, 67)= 0
pst(2, 67)= 0
pst(1, 68)= 0
pst(2, 68)= 0
pst(1, 69)= 0
pst(2, 69)= 0
pst(1, 70)= 0
pst(2, 70)= 0
pst(1, 71)= 0
pst(2, 71)= 0
pst(1, 72)= 0
pst(2, 72)= 0
pst(1, 73)= 0
pst(2, 73)= 0
pst(1, 74)= 0
pst(2, 74)= 0
pst(1, 75)= 0
pst(2, 75)= 0
pst(1, 76)= 0
pst(2, 76)= 0
pst(1, 77)= 0
pst(2, 77)= 0
pst(1, 78)= 0
pst(2, 78)= 0
pst(1, 79)= 0
pst(2, 79)= 0
pst(1, 80)= 0
pst(2, 80)= 0
pst(1, 81)= 0
pst(2, 81)= 0
pst(1, 82)= 0
pst(2, 82)= 0
pst(1, 83)= 0
pst(2, 83)= 0
pst(1, 84)= 0
pst(2, 84)= 0
pst(1, 85)= 0
pst(2, 85)= 0
pst(1, 86)= 0
pst(2, 86)= 0
pst(1, 87)= 0
pst(2, 87)= 0
pst(1, 88)= 0
pst(2, 88)= 0
pst(1, 89)= 0
pst(2, 89)= 0
pst(1, 90)= 0
pst(2, 90)= 0
pst(1, 91)= 0
pst(2, 91)= 0
pst(1, 92)= 0
pst(2, 92)= 0
pst(1, 93)= 0
pst(2, 93)= 0
pst(1, 94)= 0
pst(2, 94)= 0
pst(1, 95)= 0
pst(2, 95)= 0
pst(1, 96)= 0
pst(2, 96)= 0
pst(1, 97)= 0
pst(2, 97)= 0
pst(1, 98)= 0
pst(2, 98)= 0
pst(1, 99)= 0
pst(2, 99)= 0
pst(1,100)= 0
pst(2,100)= 0
pst(1,101)= 0
pst(2,101)= 0
pst(1,102)= 0
pst(2,102)= 0
pst(1,103)= 0
pst(2,103)= 0
pst(1,104)= 0
pst(2,104)= 0
pst(1,105)= 0
pst(2,105)= 0
pst(1,106)= 0
pst(2,106)= 0
pst(1,107)= 0
pst(2,107)= 0
pst(1,108)= 0
pst(2,108)= 0
pst(1,109)= 0
pst(2,109)= 0
pst(1,110)= 0
pst(2,110)= 0
pst(1,111)= 0
pst(2,111)= 0
pst(1,112)= 0
pst(2,112)= 0
pst(1,113)= 0
pst(2,113)= 0
pst(1,114)= 0
pst(2,114)= 0
pst(1,115)= 0
pst(2,115)= 0
pst(1,116)= 0
pst(2,116)= 0
pst(1,117)= 0
pst(2,117)= 0
pst(1,118)= 0
pst(2,118)= 0
pst(1,119)= 0
pst(2,119)= 0
pst(1,120)= 0
pst(2,120)= 0
pst(1,121)= 0
pst(2,121)= 0
pst(1,122)= 0
pst(2,122)= 0
pst(1,123)= 0
pst(2,123)= 0
pst(1,124)= 0
pst(2,124)= 0
pst(1,125)= 0
pst(2,125)= 0
pst(1,126)= 0
pst(2,126)= 0
pst(1,127)= 0
pst(2,127)= 0
pst(1,128)= 0
pst(2,128)= 0
pst(1,129)= 0
pst(2,129)= 0
pst(1,130)= 0
pst(2,130)= 0
pst(1,131)= 0
pst(2,131)= 0
pst(1,132)= 0
pst(2,132)= 0
pst(1,133)= 0
pst(2,133)= 0
pst(1,134)= 0
pst(2,134)= 0
pst(1,135)= 0
pst(2,135)= 0
pst(1,136)= 0
pst(2,136)= 0
pst(1,137)= 0
pst(2,137)= 0
pst(1,138)= 0
pst(2,138)= 0
pst(1,139)= 0
pst(2,139)= 0
pst(1,140)= 0
pst(2,140)= 0
pst(1,141)= 0
pst(2,141)= 0
pst(1,142)= 0
pst(2,142)= 0
pst(1,143)= 0
pst(2,143)= 0
pst(1,144)= 0
pst(2,144)= 0
pst(1,145)= 0
pst(2,145)= 0
pst(1,146)= 0
pst(2,146)= 0
pst(1,147)= 0
pst(2,147)= 0
pst(1,148)= 0
pst(2,148)= 0
pst(1,149)= 0
pst(2,149)= 0
pst(1,150)= 0
pst(2,150)= 0
pst(1,151)= 0
pst(2,151)= 0
pst(1,152)= 0
pst(2,152)= 0
pst(1,153)= 0
pst(2,153)= 0
pst(1,154)= 0
pst(2,154)= 0
pst(1,155)= 0
pst(2,155)= 0
pst(1,156)= 0
pst(2,156)= 0
pst(1,157)= 0
pst(2,157)= 0
pst(1,158)= 0
pst(2,158)= 0
pst(1,159)= 0
pst(2,159)= 0
pst(1,160)= 0
pst(2,160)= 0
pst(1,161)= 0
pst(2,161)= 0
pst(1,162)= 0
pst(2,162)= 0
pst(1,163)= 0
pst(2,163)= 0
pst(1,164)= 0
pst(2,164)= 0
pst(1,165)= 0
pst(2,165)= 0
pst(1,166)= 0
pst(2,166)= 0
pst(1,167)= 0
pst(2,167)= 0
pst(1,168)= 0
pst(2,168)= 0
pst(1,169)= 0
pst(2,169)= 0
pst(1,170)= 0
pst(2,170)= 0
pst(1,171)= 0
pst(2,171)= 0
pst(1,172)= 0
pst(2,172)= 0
pst(1,173)= 0
pst(2,173)= 0
pst(1,174)= 0
pst(2,174)= 0
pst(1,175)= 0
pst(2,175)= 0
pst(1,176)= 0
pst(2,176)= 0
pst(1,177)= 0
pst(2,177)= 0
pst(1,178)= 0
pst(2,178)= 0
pst(1,179)= 0
pst(2,179)= 0
pst(1,180)= 0
pst(2,180)= 0
pst(1,181)= 0
pst(2,181)= 0
pst(1,182)= 0
pst(2,182)= 0
pst(1,183)= 0
pst(2,183)= 0
pst(1,184)= 0
pst(2,184)= 0
pst(1,185)= 0
pst(2,185)= 0
pst(1,186)= 0
pst(2,186)= 0
pst(1,187)= 0
pst(2,187)= 0
pst(1,188)= 0
pst(2,188)= 0
pst(1,189)= 0
pst(2,189)= 0
pst(1,190)= 0
pst(2,190)= 0
pst(1,191)= 0
pst(2,191)= 0
pst(1,192)= 0
pst(2,192)= 0
pst(1,193)= 0
pst(2,193)= 0
pst(1,194)= 0
pst(2,194)= 0
pst(1,195)= 0
pst(2,195)= 0
pst(1,196)= 0
pst(2,196)= 0
pst(1,197)= 0
pst(2,197)= 0
pst(1,198)= 0
pst(2,198)= 0
pst(1,199)= 0
pst(2,199)= 0
pst(1,200)= 0
pst(2,200)= 0
pst(1,201)= 0
pst(2,201)= 0
pst(1,202)= 0
pst(2,202)= 0
pst(1,203)= 0
pst(2,203)= 0
pst(1,204)= 0
pst(2,204)= 0
pst(1,205)= 0
pst(2,205)= 0
pst(1,206)= 0
pst(2,206)= 0
pst(1,207)= 0
pst(2,207)= 0
pst(1,208)= 0
pst(2,208)= 0
pst(1,209)= 0
pst(2,209)= 0
pst(1,210)= 0
pst(2,210)= 0
pst(1,211)= 0
pst(2,211)= 0
pst(1,212)= 0
pst(2,212)= 0
pst(1,213)= 0
pst(2,213)= 0
pst(1,214)= 0
pst(2,214)= 0
pst(1,215)= 0
pst(2,215)= 0
pst(1,216)= 0
pst(2,216)= 0
pst(1,217)= 0
pst(2,217)= 0
pst(1,218)= 0
pst(2,218)= 0
pst(1,219)= 0
pst(2,219)= 0
pst(1,220)= 0
pst(2,220)= 0
pst(1,221)= 0
pst(2,221)= 0
pst(1,222)= 0
pst(2,222)= 0
pst(1,223)= 0
pst(2,223)= 0
pst(1,224)= 0
pst(2,224)= 0
pst(1,225)= 0
pst(2,225)= 0
pst(1,226)= 0
pst(2,226)= 0
pst(1,227)= 0
pst(2,227)= 0
pst(1,228)= 0
pst(2,228)= 0
pst(1,229)= 0
pst(2,229)= 0
pst(1,230)= 0
pst(2,230)= 0
pst(1,231)= 0
pst(2,231)= 0
pst(1,232)= 0
pst(2,232)= 0
pst(1,233)= 0
pst(2,233)= 0
pst(1,234)= 0
pst(2,234)= 0
pst(1,235)= 0
pst(2,235)= 0
pst(1,236)= 0
pst(2,236)= 0
pst(1,237)= 0
pst(2,237)= 0
pst(1,238)= 0
pst(2,238)= 0
pst(1,239)= 0
pst(2,239)= 0
pst(1,240)= 0
pst(2,240)= 0
pst(1,241)= 0
pst(2,241)= 0
pst(1,242)= 0
pst(2,242)= 0
pst(1,243)= 0
pst(2,243)= 0
pst(1,244)= 0
pst(2,244)= 0
pst(1,245)= 0
pst(2,245)= 0
pst(1,246)= 0
pst(2,246)= 0
pst(1,247)= 0
pst(2,247)= 0
pst(1,248)= 0
pst(2,248)= 0
pst(1,249)= 0
pst(2,249)= 0
pst(1,250)= 0
pst(2,250)= 0
pst(1,251)= 0
pst(2,251)= 0
pst(1,252)= 0
pst(2,252)= 0
pst(1,253)= 0
pst(2,253)= 0
pst(1,254)= 0
pst(2,254)= 0
pst(1,255)= 0
pst(2,255)= 0
pst(1,256)= 0
pst(2,256)= 0
pst(1,257)= 0
pst(2,257)= 0
pst(1,258)= 0
pst(2,258)= 0
pst(1,259)= 0
pst(2,259)= 0
pst(1,260)= 0
pst(2,260)= 0
pst(1,261)= 0
pst(2,261)= 0
pst(1,262)= 0
pst(2,262)= 0
pst(1,263)= 0
pst(2,263)= 0
pst(1,264)= 0
pst(2,264)= 0
pst(1,265)= 0
pst(2,265)= 0
pst(1,266)= 0
pst(2,266)= 0
pst(1,267)= 0
pst(2,267)= 0
pst(1,268)= 0
pst(2,268)= 0
pst(1,269)= 0
pst(2,269)= 0
pst(1,270)= 0
pst(2,270)= 0
pst(1,271)= 0
pst(2,271)= 0
pst(1,272)= 0
pst(2,272)= 0
pst(1,273)= 0
pst(2,273)= 0
pst(1,274)= 0
pst(2,274)= 0
pst(1,275)= 0
pst(2,275)= 0
pst(1,276)= 0
pst(2,276)= 0
pst(1,277)= 0
pst(2,277)= 0
pst(1,278)= 0
pst(2,278)= 0
pst(1,279)= 0
pst(2,279)= 0
pst(1,280)= 0
pst(2,280)= 0
pst(1,281)= 0
pst(2,281)= 0
pst(1,282)= 0
pst(2,282)= 0
pst(1,283)= 0
pst(2,283)= 0
pst(1,284)= 0
pst(2,284)= 0
pst(1,285)= 0
pst(2,285)= 0
pst(1,286)= 0
pst(2,286)= 0
pst(1,287)= 0
pst(2,287)= 0
pst(1,288)= 0
pst(2,288)= 0
pst(1,289)= 0
pst(2,289)= 0
pst(1,290)= 0
pst(2,290)= 0
pst(1,291)= 0
pst(2,291)= 0
pst(1,292)= 0
pst(2,292)= 0
pst(1,293)= 0
pst(2,293)= 0
pst(1,294)= 0
pst(2,294)= 0
pst(1,295)= 0
pst(2,295)= 0
pst(1,296)= 0
pst(2,296)= 0
pst(1,297)= 0
pst(2,297)= 0
pst(1,298)= 0
pst(2,298)= 0
pst(1,299)= 0
pst(2,299)= 0
pst(1,300)= 0
pst(2,300)= 0
pst(1,301)= 0
pst(2,301)= 0
pst(1,302)= 0
pst(2,302)= 0
pst(1,303)= 0
pst(2,303)= 0
pst(1,304)= 0
pst(2,304)= 0
pst(1,305)= 0
pst(2,305)= 0
pst(1,306)= 0
pst(2,306)= 0
pst(1,307)= 0
pst(2,307)= 0
pst(1,308)= 0
pst(2,308)= 0
pst(1,309)= 0
pst(2,309)= 0
pst(1,310)= 0
pst(2,310)= 0
pst(1,311)= 0
pst(2,311)= 0
pst(1,312)= 0
pst(2,312)= 0
pst(1,313)= 0
pst(2,313)= 0
pst(1,314)= 0
pst(2,314)= 0
pst(1,315)= 0
pst(2,315)= 0
pst(1,316)= 0
pst(2,316)= 0
pst(1,317)= 0
pst(2,317)= 0
pst(1,318)= 0
pst(2,318)= 0
pst(1,319)= 0
pst(2,319)= 0
pst(1,320)= 0
pst(2,320)= 0
pst(1,321)= 0
pst(2,321)= 0
pst(1,322)= 0
pst(2,322)= 0
pst(1,323)= 0
pst(2,323)= 0
pst(1,324)= 0
pst(2,324)= 0
pst(1,325)= 0
pst(2,325)= 0
pst(1,326)= 0
pst(2,326)= 0
pst(1,327)= 0
pst(2,327)= 0
pst(1,328)= 0
pst(2,328)= 0
pst(1,329)= 0
pst(2,329)= 0
pst(1,330)= 0
pst(2,330)= 0
pst(1,331)= 0
pst(2,331)= 0
pst(1,332)= 0
pst(2,332)= 0
pst(1,333)= 0
pst(2,333)= 0
pst(1,334)= 0
pst(2,334)= 0
pst(1,335)= 0
pst(2,335)= 0
pst(1,336)= 0
pst(2,336)= 0
pst(1,337)= 0
pst(2,337)= 0
pst(1,338)= 0
pst(2,338)= 0
pst(1,339)= 0
pst(2,339)= 0
pst(1,340)= 0
pst(2,340)= 0
pst(1,341)= 0
pst(2,341)= 0
pst(1,342)= 0
pst(2,342)= 0
pst(1,343)= 0
pst(2,343)= 0
pst(1,344)= 0
pst(2,344)= 0
pst(1,345)= 0
pst(2,345)= 0
pst(1,346)= 0
pst(2,346)= 0
pst(1,347)= 0
pst(2,347)= 0
pst(1,348)= 0
pst(2,348)= 0
pst(1,349)= 0
pst(2,349)= 0
pst(1,350)= 0
pst(2,350)= 0
pst(1,351)= 0
pst(2,351)= 0
pst(1,352)= 0
pst(2,352)= 0
pst(1,353)= 0
pst(2,353)= 0
pst(1,354)= 0
pst(2,354)= 0
pst(1,355)= 0
pst(2,355)= 0
pst(1,356)= 0
pst(2,356)= 0
pst(1,357)= 0
pst(2,357)= 0
pst(1,358)= 0
pst(2,358)= 0
pst(1,359)= 0
pst(2,359)= 0
pst(1,360)= 0
pst(2,360)= 0
pst(1,361)= 0
pst(2,361)= 0
pst(1,362)= 0
pst(2,362)= 0
pst(1,363)= 0
pst(2,363)= 0
pst(1,364)= 0
pst(2,364)= 0
pst(1,365)= 0
pst(2,365)= 0
pst(1,366)= 0
pst(2,366)= 0
pst(1,367)= 0
pst(2,367)= 0
pst(1,368)= 0
pst(2,368)= 0
pst(1,369)= 0
pst(2,369)= 0
pst(1,370)= 0
pst(2,370)= 0
pst(1,371)= 0
pst(2,371)= 0
pst(1,372)= 0
pst(2,372)= 0
pst(1,373)= 0
pst(2,373)= 0
pst(1,374)= 0
pst(2,374)= 0
pst(1,375)= 0
pst(2,375)= 0
pst(1,376)= 0
pst(2,376)= 0
pst(1,377)= 0
pst(2,377)= 0
pst(1,378)= 0
pst(2,378)= 0
pst(1,379)= 0
pst(2,379)= 0
pst(1,380)= 0
pst(2,380)= 0
pst(1,381)= 0
pst(2,381)= 0
pst(1,382)= 0
pst(2,382)= 0
pst(1,383)= 0
pst(2,383)= 0
pst(1,384)= 0
pst(2,384)= 0
pst(1,385)= 0
pst(2,385)= 0
pst(1,386)= 0
pst(2,386)= 0
pst(1,387)= 0
pst(2,387)= 0
pst(1,388)= 0
pst(2,388)= 0
pst(1,389)= 0
pst(2,389)= 0
pst(1,390)= 0
pst(2,390)= 0
pst(1,391)= 0
pst(2,391)= 0
pst(1,392)= 0
pst(2,392)= 0
pst(1,393)= 0
pst(2,393)= 0
pst(1,394)= 0
pst(2,394)= 0
pst(1,395)= 0
pst(2,395)= 0
pst(1,396)= 0
pst(2,396)= 0
pst(1,397)= 0
pst(2,397)= 0
pst(1,398)= 0
pst(2,398)= 0
pst(1,399)= 0
pst(2,399)= 0
pst(1,400)= 0
pst(2,400)= 0
End SUBROUTINE init_dip_planar_data
End Module dip_param

172
src/model/invariants_nh3.f90
View File

@ -1,172 +0,0 @@
Module invariants_mod
implicit none
contains
!----------------------------------------------------
subroutine invariants(a,xs,ys,xb,yb,inv)
implicit none
!include "nnparams.incl"
double precision, intent(in) :: a, xs, ys, xb, yb
double precision, intent(out) ::inv(3)
double precision:: invar(23)
!double precision:: invar(24)
complex(8) :: q1, q2
LOGICAL,PARAMETER:: debg =.FALSE.
integer :: i
! express the coordinate in complex
q1 = dcmplx(xs, ys)
q2 = dcmplx(xb, yb)
! compute the invariants
!inv( 4) = a
!inv( 4) = b**2
! INVARIANTS OF KIND II
!------------------------
invar(1) = dreal( q1 * conjg(q1) ) ! r11
invar(2) = dreal( q1 * conjg(q2) ) ! r12
invar(3) = dreal( q2 * conjg(q2) ) ! r22
invar(4) = (dimag(q1 * conjg(q2)) )**2 ! rho 12**2
!INVATIANTS OF KIND III
!------------------------
invar(5) = dreal( q1 * q1 * q1 ) ! r111
invar(6) = dreal( q1 * q1 * q2 ) ! r112
invar(7) = dreal( q1 * q2 * q2 ) ! r122
invar(8) = dreal( q2 * q2 * q2 ) ! r222
invar(9) = (dimag( q1 * q1 * q1 ))**2 ! rho111**2
invar(10) = (dimag( q1 * q1 * q2 ))**2 ! rho112 **2
invar(11) = (dimag( q1 * q2 * q2 ))**2 ! rho122**2
invar(12) = (dimag( q2 * q2 * q2 ))**2 ! rho222
! INVARIANTS OF KIND IV
!-------------------------
invar(13) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q1 ))
invar(14) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q2 ))
invar(15) = (dimag( q1 * conjg(q2)) * dimag( q1 * q2 * q2 ))
invar(16) = (dimag( q1 * conjg(q2)) * dimag( q2 * q2 * q2 ))
! INVARIANTS OF KIND V
!----------------------
invar(17) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q1 * q2 ))
invar(18) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q2 * q2 ))
invar(19) = (dimag( q1 * q1 * q1 ) * dimag( q2 * q2 * q2 ))
invar(20) = (dimag( q1 * q1 * q2 ) * dimag( q1 * q2 * q2 ))
invar(21) = (dimag( q1 * q1 * q2 ) * dimag( q2 * q2 * q2 ))
invar(22) = (dimag( q1 * q2 * q2 ) * dimag( q2 * q2 * q2 ))
invar(23) = a
! the only non zero invariant for bend pure cuts
inv(1) = invar( 1)
inv(2) = invar( 5)
inv(3) = invar( 9)
! inv(4) = invar(1)
! inv(5) = invar(5)
! inv(6) = invar(9)
!inv(17)=invar(18)
!inv(19) = invar()
!inv(2) = invar(5)
!inv(3) = invar(9)
!inv(1:22)=invar(1:22)
!inv(4:8) = invar(5:9)
!inv(9)=invar(12)
if (debg) then
write(14,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4)
write(14,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12)
write(14,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16)
write(14,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22)
write(14,*)"THE INPUT COORDINATE IN COMPLEX REPRES"
write(14,*)"---------------------------------------"
write(14,*)"xs =",dreal(q1), "ys=",dimag(q1)
endif
! modify the invariants to only consider few of them
!
!invar(13:22)=0.0d0
end subroutine invariants
! subroutine invariants_ch(a, xs, ys, xb, yb, invar)
! implicit none
! double precision, intent(in) :: a, xs, ys, xb, yb
! double precision, intent(out) :: invar(25)
!
! complex(8) :: q1, q2
! double precision :: r11, r22, r12, eta
! double precision :: I1, I2, I3, I4
!
! ! Define complex coordinates
! q1 = dcmplx(xs, ys)
! q2 = dcmplx(xb, yb)
!
! ! Quadratic invariants
! r11 = dreal(q1*conjg(q1))
! r22 = dreal(q2*conjg(q2))
! r12 = dreal(q1*conjg(q2))
! eta = dimag(q1*conjg(q2))
!
! ! Cubic imag parts (used many times)
! I1 = dimag(q1*q1*q1)
! I2 = dimag(q1*q1*q2)
! I3 = dimag(q1*q2*q2)
! I4 = dimag(q2*q2*q2)
! ! the totally symmetric invariant 25
! invar(25) = a
! ! Store invariants
! invar(1) = r11
! invar(2) = r12
! invar(3) = r22
! invar(4) = eta*eta
!
! invar(5) = dreal(q1*q1*q1)
! invar(6) = dreal(q1*q1*q2)
! invar(7) = dreal(q1*q2*q2)
! invar(8) = dreal(q2*q2*q2)
!
! invar(9) = I1*I1
! invar(10) = I2*I2
! invar(11) = I3*I3
! invar(12) = I4*I4
!
! invar(13) = eta * I1
! invar(14) = eta * I2
! invar(15) = eta * I3
! invar(16) = eta * I4
!
! ! Pure quartic invariants
! invar(17) = r11*r11 + r22*r22
! invar(18) = r11*r22 - r12*r12 - eta*eta
!
! ! Cubic-imag bilinear invariants (minimal independent set)
! invar(19) = I1*I2
! invar(20) = I1*I3
! invar(21) = I1*I4
! invar(22) = I2*I3
! invar(23) = I2*I4
!
! ! Odd cubic invariant
! invar(24) = eta**3
!
! if (debg) then
! write(14,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4)
! write(14,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12)
! write(14,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16)
! write(14,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22)
!
! write(14,*)"THE INPUT COORDINATE IN COMPLEX REPRES"
! write(14,*)"---------------------------------------"
! write(14,*)"xs =",dreal(q1), "ys=",dimag(q1)
! endif
!
! end subroutine
end module invariants_mod

87
src/model/newmodes.f
View File

@ -1,87 +0,0 @@
! Author: jnshuti
! Created: 2025-10-24 16:45:07
! Last modified: 2025-10-24 16:45:27 jnshuti
! Modes:
! Mode 1
mode( 1, 1) = 0.0000000000d0
mode( 2, 1) = 0.0000000000d0
mode( 3, 1) = 0.0000000000d0
mode( 4, 1) = 0.5773502692d0
mode( 5, 1) = 0.0000000000d0
mode( 6, 1) = 0.0000000000d0
mode( 7, 1) = -0.2886751346d0
mode( 8, 1) = 0.5000000000d0
mode( 9, 1) = 0.0000000000d0
mode(10, 1) = -0.2886751346d0
mode(11, 1) = -0.5000000000d0
mode(12, 1) = 0.0000000000d0
! Mode 2
mode( 1, 2) = 0.3121511560d0
mode( 2, 2) = 0.0000000000d0
mode( 3, 2) = 0.0000000000d0
mode( 4, 2) = -0.7756982471d0
mode( 5, 2) = 0.0000000000d0
mode( 6, 2) = 0.0000000000d0
mode( 7, 2) = -0.1939245618d0
mode( 8, 2) = 0.3358871938d0
mode( 9, 2) = 0.0000000000d0
mode(10, 2) = -0.1939245618d0
mode(11, 2) = -0.3358871938d0
mode(12, 2) = 0.0000000000d0
! Mode 3
mode( 1, 3) = 0.0000000000d0
mode( 2, 3) = 0.3121511560d0
mode( 3, 3) = 0.0000000000d0
mode( 4, 3) = 0.0000000000d0
mode( 5, 3) = 0.0000000000d0
mode( 6, 3) = 0.0000000000d0
mode( 7, 3) = 0.3358871938d0
mode( 8, 3) = -0.5817736853d0
mode( 9, 3) = 0.0000000000d0
mode(10, 3) = -0.3358871938d0
mode(11, 3) = -0.5817736853d0
mode(12, 3) = 0.0000000000d0
! Mode 4
mode( 1, 4) = 0.2830826954d0
mode( 2, 4) = 0.0000000000d0
mode( 3, 4) = 0.0000000000d0
mode( 4, 4) = 0.0759441281d0
mode( 5, 4) = 0.0000000000d0
mode( 6, 4) = 0.0000000000d0
mode( 7, 4) = -0.5655692225d0
mode( 8, 4) = -0.3703779057d0
mode( 9, 4) = 0.0000000000d0
mode(10, 4) = -0.5655692225d0
mode(11, 4) = 0.3703779057d0
mode(12, 4) = 0.0000000000d0
! Mode 5
mode( 1, 5) = 0.0000000000d0
mode( 2, 5) = 0.2830826954d0
mode( 3, 5) = 0.0000000000d0
mode( 4, 5) = 0.0000000000d0
mode( 5, 5) = -0.7794070061d0
mode( 6, 5) = 0.0000000000d0
mode( 7, 5) = -0.3703779057d0
mode( 8, 5) = -0.1378936554d0
mode( 9, 5) = 0.0000000000d0
mode(10, 5) = 0.3703779057d0
mode(11, 5) = -0.1378936554d0
mode(12, 5) = 0.0000000000d0
! Mode 6
mode( 1, 6) = 0.0000000000d0
mode( 2, 6) = 0.0000000000d0
mode( 3, 6) = 0.4213954872d0
mode( 4, 6) = 0.0000000000d0
mode( 5, 6) = 0.0000000000d0
mode( 6, 6) = -0.5235856642d0
mode( 7, 6) = 0.0000000000d0
mode( 8, 6) = 0.0000000000d0
mode( 9, 6) = -0.5235856642d0
mode(10, 6) = 0.0000000000d0
mode(11, 6) = 0.0000000000d0
mode(12, 6) = -0.5235856642d0

87
src/model/nnadia_nh3.f
View File

@ -1,87 +0,0 @@
subroutine nninit_mod()
use nn_params
implicit none
end subroutine
subroutine nnadia(coords,coeffs,adiaoutp)
USE pesmod, only: potential_nh3_low
Use diabmodel, only: diab_x
use nn_params
use nncommons
implicit none
! returns THE DIABATIC MATRIX UPPER MATRIX OF BOTH X AND Y COMPONENT OF DIPOLE
!
! coords: vector containing symmetry-adapted coordinates.
! coeffs: vector containing coeffs of diab dipole matrix
! adiaoutp: upper diab matrix of dipole
! x component come first and then y component
!
! dmat_x & dmat_y: analytic model of dipole, x and y
!
! matdim: dimension of hamiltoniam matrix
! mkeigen: integer: no eigenvectors if =0
! eigenvs: dummy storage for eigenvectors
!include 'JTmod.incl'
DOUBLE PRECISION,intent(in):: coords(maxnin)
DOUBLE PRECISION,INTENT(INOUT)::coeffs(maxnout)
DOUBLE PRECISION,INTENT(OUT):: adiaoutp(maxpout)
DOUBLE PRECISION, DIMENSION(ndiab,ndiab):: v,U,ex,temp
double precision, dimension(ndiab):: e
Double precision, allocatable:: mat(:,:)
integer i, ij, j, max_row
! variable for lapack
integer,parameter :: lwork = 1000
double precision work(lwork)
integer info
call potential_nh3_low(V,coords,coeffs)
! diagonalize the pes
u=0.0d0
allocate(mat,source=V)
call DSYEV('V','U',ndiab,mat,ndiab,e,work,lwork, info)
if (info .ne. 0) then
write(6,*)"The diagonalizing routine failed", info
stop
end if
U(:,:)= mat(:,:)
deallocate(mat)
do i =1,ndiab
adiaoutp(i)=e(i)
enddo
! Fix the phase of U matrix
! fix phas of U
do i =1,ndiab
max_row=maxloc(abs(U(:,i)),1)
if(U(max_row,i) .lt. 0) then
U(:,i) = -1.0d0*U(:,i)
endif
enddo
call diab_x(ex,coords,coeffs) ! ex is the dipole x component
! transforrm the dipole
temp = matmul(transpose(U),ex)
ex = matmul(temp,U)
ij = 4
do i =1,ndiab
do j = i,ndiab
ij = ij+1
adiaoutp(ij) = ex(i,j)
!adiaoutp(ij) = 0.0d0
enddo
enddo
END SUBROUTINE

91
src/model/nncoords_nh3.f90
View File

@ -1,91 +0,0 @@
module trans_coord
implicit none
contains
subroutine cart2mode(qq)
use invariants_mod, only: invariants
use nn_params
use nncommons
implicit none
double precision, intent(inout):: qq(maxnin)
double precision :: q(12),inv(3),qm(6)
double precision, parameter :: m_N=25526.04237518618d0
double precision, parameter :: m_H=1837.152647439619d0
double precision, parameter :: Ang2Bohr=1.889726124565d0
double precision diff(12), ref_geom(12)
double precision q_mod(12)
! transformation matrix
double precision mode(12,6)
integer i
! final internal coord vector
include 'newmodes.f'
ref_geom(1:12)=0.d0
ref_geom( 4 )=1.022871078d0
ref_geom( 7 )=-.5d0*ref_geom(4)
!ref_geom( 8 )=sqrt(1.5d0)*ref_geom(4)
ref_geom( 8)=0.5d0*sqrt(3.0d0)*ref_geom(4)
ref_geom(10 )=-.5d0*ref_geom(4)
!ref_geom(11 )=-sqrt(1.5d0)*ref_geom(4)
ref_geom(11)=-0.5d0*sqrt(3.0d0)*ref_geom(4)
!massN=25526.04237518618
!massH=1837.152647439619
!
q(1:12)=qq(1:12)
q_mod = (1.0d0/Ang2Bohr)*q
! difference vector
diff=q_mod-ref_geom
! mass weighting
diff(1:3)=diff(1:3)*sqrt(m_N)
diff(4:12)=diff(4:12)*sqrt(m_H)
qm=matmul(transpose(mode),diff)
do i =1,size(qm)
if (abs(qm(i)) .lt. 1.0d-6) qm(i)=0.0d0
enddo
!qm(3) = -qm(3)
call invariants(qm(1),qm(2),qm(3),qm(4),qm(5),inv)
qq =0.0d0
qq(1:len_in)=inv(1:len_in)
qq(len_in+1:len_in+6)=qm(1:6)
end subroutine cart2mode
end module trans_coord
!**** Define coordinate transformation applied to the input before fit.
!***
!***
!****Conventions:
!***
!*** ctrans: subroutine transforming a single point in coordinate space
subroutine trans_in(pat_in)
!use ctrans_mod, only: ctrans
use trans_coord, only: cart2mode
use nn_params
use nncommons
implicit none
!include 'nndbg.incl'
double precision pat_in(maxnin,maxpats)
!integer ntot
integer j,i,jj
jj=1
do j=1,sets
!write(62,*)"Scan ",j
do i =1,ndata(j)
call cart2mode(pat_in(:,jj))
! FIRST ELEMENT OF PAT-IN ARE USED BY NEURON NETWORK
!write(62,'(*(ES16.8))') pat_in(:,jj)
jj = jj +1
enddo
enddo
end

172
src/model/pes_model.f90
View File

@ -1,172 +0,0 @@
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

72
src/model/red_invariants.f90
View File

@ -1,72 +0,0 @@
Module invariants_mod
implicit none
contains
!----------------------------------------------------
subroutine invariants(a,xs,ys,xb,yb,b,invar)
implicit none
double precision, intent(in) :: a, xs, ys, xb, yb, b
double precision, intent(out) :: invar(24)
complex(8) :: q1, q2
LOGICAL,PARAMETER:: debg =.FALSE.
integer :: i
! express the coordinate in complex
q1 = dcmplx(xs, ys)
q2 = dcmplx(xb, yb)
! compute the invariants
invar(24) = a
invar(23) =b**2
! INVARIANTS OF KIND II
!------------------------
invar(1) = dreal( q1 * conjg(q1) ) ! r11
invar(2) = dreal( q1 * conjg(q2) ) ! r12
invar(3) = dreal( q2 * conjg(q2) ) ! r22
invar(4) = (dimag(q1 * conjg(q2)) )**2 ! rho 12**2
!INVATIANTS OF KIND III
!------------------------
invar(5) = dreal( q1 * q1 * q1 ) ! r111
invar(6) = dreal( q1 * q1 * q2 ) ! r112
invar(7) = dreal( q1 * q2 * q2 ) ! r122
invar(8) = dreal( q2 * q2 * q2 ) ! r222
invar(9) = (dimag( q1 * q1 * q1 ))**2 ! rho111**2
invar(10) = (dimag( q1 * q1 * q2 ))**2 ! rho112 **2
invar(11) = (dimag( q1 * q2 * q2 ))**2 ! rho122**2
invar(12) = (dimag( q2 * q2 * q2 ))**2 ! rho222
! INVARIANTS OF KIND IV
!-------------------------
invar(13) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q1 ))
invar(14) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q2 ))
invar(15) = (dimag( q1 * conjg(q2)) * dimag( q1 * q2 * q2 ))
invar(16) = (dimag( q1 * conjg(q2)) * dimag( q2 * q2 * q2 ))
! INVARIANTS OF KIND V
!----------------------
invar(17) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q1 * q2 ))
invar(18) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q2 * q2 ))
invar(19) = (dimag( q1 * q1 * q1 ) * dimag( q2 * q2 * q2 ))
invar(20) = (dimag( q1 * q1 * q2 ) * dimag( q1 * q2 * q2 ))
invar(21) = (dimag( q1 * q1 * q2 ) * dimag( q2 * q2 * q2 ))
invar(22) = (dimag( q1 * q2 * q2 ) * dimag( q2 * q2 * q2 ))
if (debg) then
write(*,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4)
write(*,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12)
write(*,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16)
write(*,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22)
write(*,*)"THE INPUT COORDINATE IN COMPLEX REPRES"
write(*,*)"---------------------------------------"
write(*,*)"xs =",dreal(q1), "ys=",dimag(q1)
endif
end subroutine invariants
end module invariants_mod

9
src/model/sh-scal-50.incl
View File

@ -1,9 +0,0 @@
! shift and scale for 70N structrure
! Real(dp),parameter:: scal(17) = 1.0d0 ,shift(17) = 1.0d0
Real(dp),parameter:: scal(17) = &
[0.00434,0.00239,0.00201,4.13538,9.04704,1.07732,0.79253, &
1.86682,1.87266,0.58857,0.90899,2.07755,1.81564,1.31957,1.68503,1.48391,1.79901]
Real(dp),parameter:: shift(17) = &
[-0.00016,-0.00020, 0.00008,-0.04261,13.00281,-0.44002,-0.18611, &
0.73999, 0.34618,-0.44183,-0.06997, 0.39988,-1.42376, 0.77752,-1.23819,-0.22926, 0.82597]

11
src/model/shift-scale-70N.incl
View File

@ -1,11 +0,0 @@
! shift and scale for 70N structrure
! Real(dp),parameter:: scal(17) = 1.0d0 ,shift(17) = 1.0d0
Real(dp),parameter:: scal(17) = &
[0.00401, 0.00265, 0.00335, 1.19496, 3.79876, 0.32139, 0.46267, &
1.36249, 0.90118, 0.94406, 0.29014, 0.72172, 1.25435, 2.09880, &
0.93696, 1.69131, 0.75050 ]
Real(dp),parameter:: shift(17) = &
[-0.00018,-0.00018,-0.00002,-0.03463,-0.74584, 0.09806, 0.16350, &
-0.96733,-0.62850, 0.57182,-0.00110, 0.10584, 0.29949,-0.82418, &
0.57450, 0.04471, 0.28142 ]

235
src/model/write_error.f90
View File

@ -1,235 +0,0 @@
module print_error
use iso_fortran_env, only: idp => int32, dp => real64
use nn_params
use nncommons
implicit none
contains
!=======================================================================
! Compute and print error summary (RMS) for PES energies and dipole moments
! Separately for total, weighted, energies (1:4), dipoles (5:inp_out)
!=======================================================================
subroutine print_ErrorSummary(par, pat_in, pat_out,ref_in,ref_out, &
wterr, typop, laystr, weistr, nlay, &
npat, nref)
real(dp), intent(in) :: pat_in(maxnin, maxpats)
real(dp), intent(in) :: pat_out(maxpout, maxpats)
real(dp), intent(in) :: par(wbcap, maxset)
real(dp), intent(in) :: wterr(maxpout, maxpats)
real(dp), intent(in) :: ref_in(maxnin,maxpats), ref_out(maxpout,maxpats)
integer(idp), intent(in) :: laystr(3, maxlay)
integer(idp), intent(in) :: typop(maxtypes, maxlay)
integer(idp), intent(in) :: weistr(2, maxlay, 2)
integer(idp), intent(in) :: nlay, npat, nref
! Local variables
integer(idp) :: pt, k,i
integer(idp) :: ntot
integer(idp),dimension(maxpout) :: cnt_tot,cnt_fit,cnt_val
real(dp), dimension(maxpout) :: total_rms_fit, total_rms_val,wrms_fit,wrms_val
real(dp) :: nnoutp(maxnout)
real(dp) :: ymod(maxpout, maxpats)
real(dp) :: wsum_fit
real(dp) :: rms_total, rms_weight
real(dp) :: en_rms, dip_rms
Real(dp) :: rms_total_fit, rms_w_fit
Real(dp) :: rms_total_val, rms_w_val
Real(dp) :: en_rms_fit,en_rms_val,dip_rms_fit,dip_rms_val
integer(idp), parameter :: std_out = 6
! Total number of points (fitting + validation/reference)
ntot = npat + nref
! keep in mind that npat is the total number of point fit + ref
!write(6,*)"Npat in the error summary", npat,nref
write(11,'(*(ES16.7))')par(:,1)
! Initialize accumulators
cnt_fit = 0
cnt_tot = 0
cnt_val = 0
wsum_fit =0.0_dp
total_rms_fit(:) = 0.0_dp
wrms_fit(:) = 0.0_dp
total_rms_val (:) = 0.0_dp
wrms_val(:) = 0.0_dp
!write(11,'(*(ES18.6))')ref_in(1:3,:)
! Generate model predictions and accumulate squared errors
pt = 0
do pt =1, npat
if (pt > maxpats) then
write(std_out, *) "ERROR: pt exceeds maxpats!"
stop
end if
! Get network output
call neunet(pat_in(:, pt), nnoutp, par,typop, laystr, weistr, nlay)
! Transform to adiabatic representation (energies + dipoles)
call nnadia(pat_in(:, pt), nnoutp, ymod(:, pt))
! Accumulate errors only for components with positive weight
do k = 1, inp_out
if (wterr(k, pt) > 0.0_dp) then
!cnt_tot(k) = cnt_tot(k) + 1
! Fitting / training set
cnt_fit(k) = cnt_fit(k) + 1
total_rms_fit(k) = total_rms_fit(k) + (ymod(k,pt) - pat_out(k,pt))**2
wrms_fit(k) = wrms_fit(k) + ((ymod(k,pt) - pat_out(k,pt))**2 &
* (wterr(k,pt)))
wsum_fit = wsum_fit + wterr(k,pt)
end if
end do
end do
! validation set
do i =1, nref
!if (pt > maxpats) stop "Error pt exceeds maxpats"
pt = npat + i
call neunet(ref_in(:,i),nnoutp,par,typop,laystr,weistr,nlay)
! call the model
call nnadia(ref_in(:,i),nnoutp,ymod(:,pt))
! calculate the errors
do k=1,inp_out
if (wterr(k,pt) > 0.0_dp) then
cnt_val (k) = cnt_val(k) + 1
total_rms_val(k) = total_rms_val(k) + (ymod(k,pt) - ref_out(k,i))**2
wrms_val(k) = wrms_val(k) + ((ymod(k,pt) - ref_out(k,i))**2 * &
(wterr(k,pt)))
!write(6,*)" IN here wter"
endif
enddo
!write(11,'(*(ES18.6))')wterr(:,pt)
enddo
! Safety check: did we process the expected number of points?
! Compute RMS values
! TOTAL FIT +REF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
rms_total = 0.0_dp
if ((sum(cnt_fit(1:inp_out)) + sum(cnt_val(1:inp_out))) > 0) then
rms_total = dsqrt( sum(total_rms_fit + total_rms_val) / &
real((sum(cnt_fit(1:inp_out))+sum(cnt_val(1:inp_out))), dp) )
end if
rms_weight = 0.0_dp
if (sum(wterr(1:inp_out, 1:npat+nref)) > 1.0e-12_dp) then ! avoid div-by-zero
rms_weight = dsqrt( sum(wrms_fit) + sum(wrms_val)) / sum(wterr(1:inp_out, 1:npat+nref))
!rms_weight = dsqrt(wt_rms /weighted_wt )
end if
en_rms = 0.0_dp
if ((sum(cnt_fit(1:4)) + sum(cnt_val(1:4))) > 0) then
en_rms = sqrt( ( sum(total_rms_fit(1:4)) + sum(total_rms_val(1:4)) ) / &
real(sum(cnt_fit(1:4)) + sum(cnt_val(1:4)), dp) )
end if
dip_rms = 0.0_dp
if ((sum(cnt_fit(5:inp_out)) + sum(cnt_val(5:inp_out)))> 0) then
dip_rms = sqrt( (sum(total_rms_fit(5:inp_out)) + sum(total_rms_val(5:inp_out))) / &
real((sum(cnt_fit(5:inp_out))+sum(cnt_val(5:inp_out))), dp) )
end if
rms_total_fit = 0.0_dp
if (sum(cnt_fit(1:inp_out)) > 0) then
rms_total_fit = dsqrt( sum(total_rms_fit(1:inp_out)) / real(sum(cnt_fit(1:inp_out)), dp) )
endif
rms_w_fit = 0.0_dp
if (wsum_fit > 1.0e-12_dp) then
rms_w_fit = sqrt( sum(wrms_fit(1:inp_out)) / wsum_fit )
endif
! validation
rms_total_val = 0.0_dp
if (sum(cnt_val(1:inp_out)) > 0) then
rms_total_val = sqrt( sum(total_rms_val(1:inp_out)) / real(sum(cnt_val(1:inp_out)), dp) )
endif
rms_w_val = 0.0_dp
if (sum(wterr(1:inp_out,npat+1:npat+nref)) > 1.0e-12_dp) then
rms_w_val = sqrt( sum(wrms_val(1:inp_out)) / sum(wterr(1:inp_out,npat+1:npat+nref)) )
endif
! compute each quantity error
en_rms_fit = merge( sqrt(sum(total_rms_fit(1:4))/real(sum(cnt_fit(1:4)),dp)),&
0.0_dp, sum(cnt_fit(1:4))>0 )
en_rms_val = merge( sqrt(sum(total_rms_val(1:4))/real(sum(cnt_val(1:4)),dp)), &
0.0_dp, sum(cnt_val(1:4))>0 )
dip_rms_fit = merge( sqrt(sum(total_rms_fit(5:inp_out))/real(sum(cnt_fit(5:inp_out)),dp)), &
0.0_dp, sum(cnt_fit(5:inp_out))>0 )
dip_rms_val = merge( sqrt(sum(total_rms_val(5:inp_out))/real(sum(cnt_val(5:inp_out)),dp)), &
0.0_dp, sum(cnt_val(5:inp_out))>0 )
! Output
write(std_out, '(A)') " "
write(std_out, '(A)') "USER DEFINED ERROR METRIC"
write(std_out, '(A)') "######################################"
!write(std_out, '(A,ES18.8,A)') &
! "Total RMS (all fitted components, unweighted) : ", rms_total * unit_con, &
! " [ha]"
!write(std_out, '(A,ES18.8,A)') &
! "Total Weighted RMS (fit + Ref) : ", rms_weight * unit_con
!write(std_out, '(A,ES18.8,A)') &
! "Energy RMS : ", en_rms * unit_con, " [ha] "
!write(std_out, '(A,ES18.8,A)') &
! "Dipole RMS : ", dip_rms * unit_con, " [ha] "
! compute the metric for fit data
write(std_out, '(A)') "FITTING / TRAINING SET"
write(std_out, '(A)') "#############################"
write(std_out,'(A,ES14.6,A)') " Total RMS (unweighted) = ", rms_total_fit * unit_con
write(std_out,'(A,ES14.6,A)') " Weighted RMS = ", rms_w_fit * unit_con
write(std_out,'(A,ES14.6,A)') " Energy RMS = ", en_rms_fit * unit_con
write(std_out,'(A,ES14.6,A)') " Dipole RMS = ", dip_rms_fit * unit_con
write(std_out, '(A)') " "
write(std_out, '(A)') "VALIDATION / REFERENCE SET"
write(std_out, '(A)') "##################################"
write(std_out,'(A,ES14.6,A)') " Total RMS (unweighted) = ", rms_total_val * unit_con
write(std_out,'(A,ES14.6,A)') " Weighted RMS = ", rms_w_val * unit_con
write(std_out,'(A,ES14.6,A)') " Energy RMS = ", en_rms_val * unit_con
write(std_out,'(A,ES14.6,A)') " Dipole RMS = ", dip_rms_val * unit_con
write(std_out, '(A)') " "
write(std_out, '(A)') "COMBINED (FIT + VALIDATION)"
write(std_out, '(A)') "------------------------------"
write(std_out,'(A,ES14.6,A)') " Total RMS (unweighted) = ", rms_total * unit_con
write(std_out,'(A,ES14.6,A)') " Weighted RMS = ", rms_weight * unit_con
write(std_out,'(A,ES14.6,A)') " Energy RMS = ", en_rms * unit_con
write(std_out,'(A,ES14.6,A)') " Dipole RMS = ", dip_rms * unit_con
write(std_out,'(A)')" END OF ERROR SUMMARY"
!write(std_out, '(A)') repeat("-", 70)
! Optional: add more detailed per-component output if desired
! do k = 1, inp_out
! if (cnt(k) > 0) then
! write(std_out,'(A,I0,A,ES12.5)') "Output ",k," RMS = ", &
! sqrt(total_rms(k)/real(cnt(k),dp)) * unit_con
! end if
! end do
end subroutine print_ErrorSummary
end module print_error