jnshuti
/
Fitting-L-matrix
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Update the keys

This commit is contained in:
jean paul nshuti 2025-10-06 13:48:56 +02:00
parent fc9159bc32
commit 05265f77a3
2 changed files with 102 additions and 130 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
src/model/keys.f90
View File

@ -10,7 +10,7 @@ contains
character(len=1) prefix(4) character(len=1) prefix(4)
parameter (prefix=['N','P','A','S']) parameter (prefix=['N','P','A','S'])
!character (len=20) key(4,25) !character (len=20) key(4,25)
integer,parameter:: np=34 integer,parameter:: np=32
character(len=16) parname(np) character(len=16) parname(np)
integer i,j integer i,j
! Defining keys for potential ! Defining keys for potential
@ -24,59 +24,57 @@ contains
parname(4)='LXYVA2O2' parname(4)='LXYVA2O2'
parname(5)='LXYVE1O2' parname(5)='LXYVE1O2'
parname(6)='LXYVE2O2' parname(6)='LXYVE2O2'
parname(7)='LXYVA2O3'
parname(8)='LXYVE1O3'
parname(9)='LXYVE2O3'
! W & Z of E1 ! W & Z of E1
parname(10)='LXYWZE1O0' parname(7)='LXYWZE1O0'
parname(11)='LXYWZE1O1' parname(8)='LXYWZE1O1'
parname(12)='LXYWZE1O2' parname(9)='LXYWZE1O2'
parname(13)='LXYWZE1O3' parname(10)='LXYWZE2O0'
parname(14)='LXYWZE2O0' parname(11)='LXYWZE2O1'
parname(15)='LXYWZE2O1' parname(12)='LXYWZE2O2'
parname(16)='LXYWZE2O2'
parname(17)='LXYWZE2O3'
! WW and Z Pseudo between E1 and E2 ! WW and Z Pseudo between E1 and E2
! p STANDS FOR PSEUDO JAHN-TELLER ! p STANDS FOR PSEUDO JAHN-TELLER
parname(18)='LXYPE1E2O0' parname(13)='LXYPE1E2O0'
parname(19)='LXYPE1E2O1' parname(14)='LXYPE1E2O1'
parname(20)='LXYPE1E2O2' parname(15)='LXYPE1E2O2'
! no order 3 ! no order 3
! PSEUDO A2 & E1 ! PSEUDO A2 & E1
parname(21)='LXYPA2E1O0' parname(16)='LXYPA2E1O0'
parname(22)='LXYPA2E1O1' parname(17)='LXYPA2E1O1'
parname(23)='LXYPA2E1O2' parname(18)='LXYPA2E1O2'
! Pseudo JAHN-TELLER BETWEEN A2 AND E1 ! Pseudo JAHN-TELLER BETWEEN A2 AND E1
parname(24)='LXYPA2E2O0' parname(19)='LXYPA2E2O0'
parname(25)='LXYPA2E2O1' parname(20)='LXYPA2E2O1'
parname(26)='LXYPA2E2O2' parname(21)='LXYPA2E2O2'
parname(27)='LXYPA2E2O3'
! keys for lz ! keys for lz
parname(28)='LZWZE1O1' parname(22)='LZWZE1O1'
parname(29)='LZWZE1O2' parname(23)='LZWZE1O2'
parname(30)='LZWZE2O1' parname(24)='LZWZE2O1'
parname(31)='LZWZE2O2' parname(25)='LZWZE2O2'
parname(33)='LZPE1E2O0' parname(26)='LZPE1E2O0'
parname(34)='LZPE1E2O1' parname(27)='LZPE1E2O1'
parname(35)='LZPE1E2O2' parname(28)='LZPE1E2O2'
parname(36)='LZPA2E1O1' parname(29)='LZPA2E1O1'
parname(37)='LZPA2E2O2' parname(30)='LZPA2E2O2'
parname(39)='LZPA2E2O1' parname(31)='LZPA2E2O1'
parname(34)='TYPE_CAL'! TYPE OF THE CALCULATION WHETHER IT IS THE TRACE OR SOMETHING ELSE parname(32)='TYPE_CAL'! TYPE OF THE CALCULATION WHETHER IT IS THE TRACE OR SOMETHING ELSE
do i=1,np do i=1,np
do j=1,4 do j=1,4
key(j, i)=prefix(j)//trim(parname(i))//':' ! first 86 keys are the potential keys key(j, i)=prefix(j)//trim(parname(i))//':'
write(8,*) key(j,i)
enddo enddo
write(8,*) ''
enddo enddo

164
src/model/model.f90
View File

@ -6,6 +6,8 @@
module diab_mod: module diab_mod:
use accuracy_constants, only: dp, idp use accuracy_constants, only: dp, idp
use dim_parameter, only: ndiab, nstat, ntot,npar,qn,pst use dim_parameter, only: ndiab, nstat, ntot,npar,qn,pst
private
public :: Lx_diab, Ly_diab, Lz_diab
implicit none implicit none
contains contains
@ -15,9 +17,11 @@ module diab_mod:
real(dp),dimension(:),intent(in):: q,t real(dp),dimension(:),intent(in):: q,t
real(dp),dimension(npar),intent(in):: p real(dp),dimension(npar),intent(in):: p
real(dp):: xs,ys,xb,yb,a,b real(dp):: xs,ys,xb,yb,a,b
real(dp):: v3_vec(8) real(dp):: v3_vec(8), v2(6)
integer(idp):: i,j,id integer(idp):: i,j,id
! check the dimension of the matrix ! check the dimension of the matrix
if (size(E,1) .ne. ndiab) then if (size(E,1) .ne. ndiab) then
write(*,*) " Error in Lx_diab: wrong dimension of L matrix ", size(E,1) write(*,*) " Error in Lx_diab: wrong dimension of L matrix ", size(E,1)
@ -26,7 +30,24 @@ module diab_mod:
! rewrite the coordinate array q into symmetry adapted coordinates ! rewrite the coordinate array q into symmetry adapted coordinates
call rewrite_coord(q,a,xs,ys,xb,yb,b,1) call rewrite_coord(q,a,xs,ys,xb,yb,b,1)
v2(1)=xs**2-ys**2
v2(2)=xb**2-yb**2
v2(3)=xs*xb-ys*yb
v2(4)=2*xs*ys
v2(5)=2*xb*yb
v2(6)=xs*yb+xb*ys
v3( 1) = xs*(xs**2-3*ys**2)
v3( 2) = xb*(xb**2-3*yb**2)
v3( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys
v3( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb
v3( 5) = ys*(3*xs**2-ys**2)
v3( 6) = yb*(3*xb**2-yb**2)
v3( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys
v3( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb
e = 0.0_dp
id = 1
e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 param e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 param
id=id+1 ! 2 id=id+1 ! 2
e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 p e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 p
@ -50,43 +71,27 @@ module diab_mod:
+p(pst(1,id)+2)*(xs*xb-ys*yb) +p(pst(1,id)+2)*(xs*xb-ys*yb)
e(5,5)=e(5,5)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & e(5,5)=e(5,5)+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)+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 ! 2 param
id=id+1 ! 8
e(2,2)=e(2,2)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb *sb ! 2 p
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 ! 2 p
e(5,5)=e(5,5)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb
! W and Z term of E1 ! W and Z term of E1
! order 0 ! order 0
id=id+1 ! 10 id=id+1 ! 7
e(2,2)=e(2,2)+p(pst(1,id)) e(2,2)=e(2,2)+p(pst(1,id))
e(3,3)=e(3,3)-p(pst(1,id)) e(3,3)=e(3,3)-p(pst(1,id))
!e(2,3)=e(2,3) !e(2,3)=e(2,3)
! order 1 ! order 1
id=id+1 ! 11 ! 2 param id=id+1 ! 8 ! 2 param
e(2,2)=e(2,2)+ p(pst(1,id))*xs+p(pst(1,id)+1)*xb 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(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 e(2,3)=e(2,3)- p(pst(1,id))*ys -p(pst(1,id)+1)*yb
! order 2 ! order 2
id=id+1 ! 12 ! 3p id=id+1 ! 9 ! 3p
do i=1,3 do i=1,3
e(2,2)=e(2,2)+p(pst(1,id)+(i-1))*v2(i) e(2,2)=e(2,2)+p(pst(1,id)+(i-1))*v2(i)
e(3,3)=e(3,3)-p(pst(1,id)+(i-1))*v2(i) e(3,3)=e(3,3)-p(pst(1,id)+(i-1))*v2(i)
e(2,3)=e(2,3)+ p(pst(1,id)+(i-1))*v2(i+3) e(2,3)=e(2,3)+ p(pst(1,id)+(i-1))*v2(i+3)
enddo enddo
! order 3 ! order 3
id=id+1 ! 13 ! 8 param
do i=1,4
e(2,2)=e(2,2)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i)
e(3,3)=e(3,3)-(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i)
e(2,3)=e(2,3)+(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i+4)
enddo
! try the testing of higher order terms ! try the testing of higher order terms
!e(2,3)=e(2,3)- p(pst(1,id))*ys*ss +p(pst(1,id)+1)*ss*2*xs*ys !e(2,3)=e(2,3)- p(pst(1,id))*ys*ss +p(pst(1,id)+1)*ss*2*xs*ys
@ -97,32 +102,24 @@ module diab_mod:
! W and Z for E2 ! W and Z for E2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
id=id+1 ! 14 id=id+1 ! 10
e(4,4)=e(4,4)+p(pst(1,id)) e(4,4)=e(4,4)+p(pst(1,id))
e(5,5)=e(5,5)-p(pst(1,id)) e(5,5)=e(5,5)-p(pst(1,id))
e(4,5)=e(4,5) e(4,5)=e(4,5)
! order 1 ! order 1
id=id+1 ! 2 param 15 id=id+1 !112 param 15
e(4,4)=e(4,4)+ p(pst(1,id))*xs+p(pst(1,id)+1)*xb e(4,4)=e(4,4)+ p(pst(1,id))*xs+p(pst(1,id)+1)*xb
e(5,5)=e(5,5)- (p(pst(1,id))*xs+p(pst(1,id)+1)*xb) e(5,5)=e(5,5)- (p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
e(4,5)=e(4,5)- p(pst(1,id))*ys-p(pst(1,id)+1)*yb e(4,5)=e(4,5)- p(pst(1,id))*ys-p(pst(1,id)+1)*yb
! order 2 ! order 2
id=id+1 ! 16 ! 3p id=id+1 ! 12 ! 3p
do i=1,3 do i=1,3
e(4,4)=e(4,4)+p(pst(1,id)+(i-1))*v2(i) e(4,4)=e(4,4)+p(pst(1,id)+(i-1))*v2(i)
e(5,5)=e(5,5)-p(pst(1,id)+(i-1))*v2(i) e(5,5)=e(5,5)-p(pst(1,id)+(i-1))*v2(i)
e(4,5)=e(4,5)+ p(pst(1,id)+(i-1))*v2(i+3) e(4,5)=e(4,5)+ p(pst(1,id)+(i-1))*v2(i+3)
enddo enddo
! order 3
id=id+1 ! 17 ! 8 param
do i=1,4
e(4,4)=e(4,4)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i)
e(5,5)=e(5,5)-(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i)
e(4,5)=e(4,5)+(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i+4)
enddo
! make the dipole E = b* E ! make the dipole E = b* E
e = b * e e = b * e
@ -131,20 +128,20 @@ module diab_mod:
! WW and ZZ ! WW and ZZ
id =id+1 ! 18 id =id+1 ! 13
e(2,4)=e(2,4)+p(pst(1,id))*b e(2,4)=e(2,4)+p(pst(1,id))*b
e(3,5)=e(3,5)-p(pst(1,id))*b e(3,5)=e(3,5)-p(pst(1,id))*b
! ORDER 1 ! ORDER 1
id=id+1 ! 19 ! 6 parama id=id+1 ! 14 ! 6 parama
e(2,4)=e(2,4)+b*((p(pst(1,id))+p(pst(1,id)+1)+p(pst(1,id)+2))*xs+(p(pst(1,id)+3)+p(pst(1,id)+4)+p(pst(1,id)+5))*xb) e(2,4)=e(2,4)+b*((p(pst(1,id))+p(pst(1,id)+1)+p(pst(1,id)+2))*xs+(p(pst(1,id)+3)+p(pst(1,id)+4)+p(pst(1,id)+5))*xb)
e(3,5)=e(3,5)+b*((p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb) e(3,5)=e(3,5)+b*((p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb)
e(2,5)=e(2,5)+b*((p(pst(1,id))-p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(p(pst(1,id)+3)-p(pst(1,id)+4)-p(pst(1,id)+5))*yb) e(2,5)=e(2,5)+b*((p(pst(1,id))-p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(p(pst(1,id)+3)-p(pst(1,id)+4)-p(pst(1,id)+5))*yb)
e(3,4)=e(3,4)+b*((-p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*yb) e(3,4)=e(3,4)+b*((-p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*yb)
! order 2 ! order 2
id=id+1 ! 20 id=id+1 ! 15
do i=1,3 ! 9 param do i=1,3 ! param
e(2,4)=e(2,4)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i) e(2,4)=e(2,4)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i)
e(3,5)=e(3,5)+b*(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i) e(3,5)=e(3,5)+b*(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i)
e(2,5)=e(2,5)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3) e(2,5)=e(2,5)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)
@ -155,17 +152,17 @@ module diab_mod:
! ################################################## ! ##################################################
!################################################### !###################################################
! order 0 ! order 0
id=id+1 ! 1 param ! 21 id=id+1 ! 1 param ! 16
e(1,3)=e(1,3)+b*(p(pst(1,id))) e(1,3)=e(1,3)+b*(p(pst(1,id)))
! order 1 ! order 1
id = id +1 ! 22 id = id +1 ! 17
e(1,2)=e(1,2)-b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) 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) e(1,3)=e(1,3)+b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
! order 2 ! order 2
id=id+1 ! 23 id=id+1 ! 18
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,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) & e(1,3)=e(1,3)+b*(p(pst(1,id))*(xs**2-ys**2) + p(pst(1,id)+1)*(xb**2-yb**2) &
@ -176,28 +173,21 @@ module diab_mod:
!########################################################################################################## !##########################################################################################################
! order 0 ! order 0
id =id+1 !24 id =id+1 !19
e(1,5)=e(1,5)+p(pst(1,id)) e(1,5)=e(1,5)+p(pst(1,id))
! order 1 ! order 1
id = id +1 ! 25 id = id +1 ! 20
e(1,4)=e(1,4)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) e(1,4)=e(1,4)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
e(1,5)=e(1,5)+(p(pst(1,id))*xs + p(pst(1,id)+1)*xb) e(1,5)=e(1,5)+(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
! order 2 ! order 2
id=id+1 ! 26 id=id+1 ! 21
e(1,4)=e(1,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(1,4)=e(1,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(1,5)=e(1,5)+p(pst(1,id))*(xs**2-ys**2) + p(pst(1,id)+1)*(xb**2-yb**2) & e(1,5)=e(1,5)+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)+2)*(xs*xb-ys*yb)
! order 3
id=id+1 ! 27 ! 8 param
do i=1,4
e(1,4)=e(1,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3(i+4)
e(1,5)=e(1,5)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i)
enddo
e(1,4:5) = b* e(1,4:5) e(1,4:5) = b* e(1,4:5)
@ -222,9 +212,12 @@ module diab_mod:
! rewrite the coordinate array q into symmetry adapted coordinates ! rewrite the coordinate array q into symmetry adapted coordinates
call rewrite_coord(q,a,xs,ys,xb,yb,b,1) call rewrite_coord(q,a,xs,ys,xb,yb,b,1)
e=0.0d0 v2(1)=xs**2-ys**2
ss=xs**2+ys**2 ! totaly symmetric term v2(2)=xb**2-yb**2
sb=xb**2+yb**2 v2(3)=xs*xb-ys*yb
v2(4)=2*xs*ys
v2(5)=2*xb*yb
v2(6)=xs*yb+xb*ys
v3( 1) = xs*(xs**2-3*ys**2) v3( 1) = xs*(xs**2-3*ys**2)
v3( 2) = xb*(xb**2-3*yb**2) v3( 2) = xb*(xb**2-3*yb**2)
@ -233,9 +226,10 @@ module diab_mod:
v3( 5) = ys*(3*xs**2-ys**2) v3( 5) = ys*(3*xs**2-ys**2)
v3( 6) = yb*(3*xb**2-yb**2) v3( 6) = yb*(3*xb**2-yb**2)
v3( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys v3( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys
v3( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb v3( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb
e=0.0d0 e = 0.0_dp
! V-term ! V-term
id=1 ! 1 id=1 ! 1
@ -256,87 +250,65 @@ module diab_mod:
id=id+1 ! 6 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)) 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))
e(5,5)=e(5,5)-(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(5,5)=e(5,5)-(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 ! 2
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 ! 3
e(4,4)=e(4,4)+p(pst(1,id))*ys*ss + p(pst(1,id)+1)*yb*sb
e(5,5)=e(5,5)+p(pst(1,id))*ys*ss + p(pst(1,id)+1)*yb*sb
! W and Z of E1 ! W and Z of E1
! order 0 ! order 0
id=id+1 ! 10 id=id+1 ! 7
e(2,3)=e(2,3)+p(pst(1,id)) e(2,3)=e(2,3)+p(pst(1,id))
! order 1 ! order 1
id=id+1 ! id=id+1 ! 8
e(2,2)=e(2,2)-p(pst(1,id))*ys -p(pst(1,id)+1)*yb 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(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 e(2,3)=e(2,3)-p(pst(1,id))*xs -p(pst(1,id)+1)*xb
! order 2 ! order 2
id=id+1 ! 12 id=id+1 ! 9
do i=1,3 do i=1,3
e(2,2)=e(2,2)+p(pst(1,id)+(i-1))*v2(i+3) e(2,2)=e(2,2)+p(pst(1,id)+(i-1))*v2(i+3)
e(3,3)=e(3,3)-p(pst(1,id)+(i-1))*v2(i+3) e(3,3)=e(3,3)-p(pst(1,id)+(i-1))*v2(i+3)
e(2,3)=e(2,3)-p(pst(1,id)+(i-1))*v2(i) e(2,3)=e(2,3)-p(pst(1,id)+(i-1))*v2(i)
enddo enddo
id=id+1 ! 8
do i=1,4
e(2,2)=e(2,2)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3(i+4)
e(3,3)=e(3,3)-(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3(i+4)
e(2,3)=e(2,3)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i)
enddo
!! W and Z of E2 !! W and Z of E2
!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!
! order 0 ! order 0
id=id+1 ! 14 id=id+1 ! 10
e(4,5)=e(4,5)+p(pst(1,id)) e(4,5)=e(4,5)+p(pst(1,id))
! order 1 ! order 1
id=id+1 ! 15 id=id+1 ! 11
e(4,4)=e(4,4)-p(pst(1,id))*ys -p(pst(1,id)+1)*yb e(4,4)=e(4,4)-p(pst(1,id))*ys -p(pst(1,id)+1)*yb
e(5,5)=e(5,5)+p(pst(1,id))*ys+ p(pst(1,id)+1)*yb e(5,5)=e(5,5)+p(pst(1,id))*ys+ p(pst(1,id)+1)*yb
e(4,5)=e(4,5)-p(pst(1,id))*xs -p(pst(1,id)+1)*xb e(4,5)=e(4,5)-p(pst(1,id))*xs -p(pst(1,id)+1)*xb
! order 2 ! order 2
id=id+1 ! 16 id=id+1 ! 12
do i=1,3 do i=1,3
e(4,4)=e(4,4)+p(pst(1,id)+(i-1))*v2(i+3) e(4,4)=e(4,4)+p(pst(1,id)+(i-1))*v2(i+3)
e(5,5)=e(5,5)-p(pst(1,id)+(i-1))*v2(i+3) e(5,5)=e(5,5)-p(pst(1,id)+(i-1))*v2(i+3)
e(4,5)=e(4,5)-p(pst(1,id)+(i-1))*v2(i) e(4,5)=e(4,5)-p(pst(1,id)+(i-1))*v2(i)
enddo enddo
id=id+1 ! 17
do i=1,4
e(4,4)=e(4,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3(i+4)
e(5,5)=e(5,5)-(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3(i+4)
e(4,5)=e(4,5)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i)
enddo
! PSEUDO JAHN-TELLER E1 AND E2 ! PSEUDO JAHN-TELLER E1 AND E2
e = b* e e = b* e
!ORDER 0 !ORDER 0
id=id+1 ! 18 id=id+1 ! 13
e(2,5)=e(2,5)+p(pst(1,id)) e(2,5)=e(2,5)+p(pst(1,id))
e(3,4)=e(3,4)+p(pst(1,id)) e(3,4)=e(3,4)+p(pst(1,id))
! order 1 ! order 1
id=id+1 id=id+1 ! 14
e(2,4)=e(2,4)+((p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*yb) e(2,4)=e(2,4)+((p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*yb)
e(3,5)=e(3,5)+((p(pst(1,id))+p(pst(1,id)+1)+p(pst(1,id)+2))*ys+(p(pst(1,id)+3)+p(pst(1,id)+4)+p(pst(1,id)+5))*yb) e(3,5)=e(3,5)+((p(pst(1,id))+p(pst(1,id)+1)+p(pst(1,id)+2))*ys+(p(pst(1,id)+3)+p(pst(1,id)+4)+p(pst(1,id)+5))*yb)
e(2,5)=e(2,5)+((-p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb) e(2,5)=e(2,5)+((-p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb)
e(3,4)=e(3,4)+((p(pst(1,id))-p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(+p(pst(1,id)+3)-p(pst(1,id)+4)-p(pst(1,id)+5))*xb) e(3,4)=e(3,4)+((p(pst(1,id))-p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(+p(pst(1,id)+3)-p(pst(1,id)+4)-p(pst(1,id)+5))*xb)
! order 2 ! order 2
id=id+1 id=id+1 ! 15
e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3) e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)
e(3,5)=e(3,5)+(-p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3) e(3,5)=e(3,5)+(-p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)
@ -350,16 +322,16 @@ module diab_mod:
! ##################### ! #####################
! order 0 ! order 0
id=id+1 id=id+1 ! 16
e(1,2)=e(1,2)+(p(pst(1,id))) e(1,2)=e(1,2)+(p(pst(1,id)))
! order 1 ! order 1
id=id+1 id=id+1 ! 17
e(1,2)=e(1,2)-(p(pst(1,id))*xs + p(pst(1,id)+1)*xb) e(1,2)=e(1,2)-(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
e(1,3)=e(1,3)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) e(1,3)=e(1,3)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
! order 2 ! order 2
id=id+1 id=id+1 !18
e(1,2)=e(1,2)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & e(1,2)=e(1,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)+2)*(xs*xb-ys*yb))
e(1,3)=e(1,3)+(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+ & e(1,3)=e(1,3)+(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+ &
@ -370,16 +342,16 @@ module diab_mod:
!############################################################################### !###############################################################################
! order 0 ! order 0
id = id+1 id = id+1 !19
e(1,4)=e(1,4)+p(pst(1,id)) e(1,4)=e(1,4)+p(pst(1,id))
! order 1 ! order 1
id=id+1 id=id+1 ! 20
e(1,4)=e(1,4)-(p(pst(1,id))*xs + p(pst(1,id)+1)*xb) e(1,4)=e(1,4)-(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
e(1,5)=e(1,5)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) e(1,5)=e(1,5)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
! order 2 ! order 2
id=id+1 id=id+1 ! 21
e(1,4)=e(1,4)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & e(1,4)=e(1,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)+2)*(xs*xb-ys*yb))
e(1,5)=e(1,5)+(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+ & e(1,5)=e(1,5)+(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+ &
@ -408,9 +380,11 @@ module diab_mod:
stop stop
endif endif
call rewrite_coord(q,a,xs,xb,yb,b,1) call rewrite_coord(q,a,xs,xb,yb,b,1)
e = 0.0_dp e = 0.0_dp
! id for lz ! id for lz
id = ! has to be id = 22 ! has to be
! the diagonal terms ! the diagonal terms
! the v-term is 0th order and 3rd order. ! the v-term is 0th order and 3rd order.
@ -418,7 +392,7 @@ module diab_mod:
! w and z of E'' ! w and z of E''
! order 1 ! order 1
id =id +1 id = id
e(2,2) = e(2,2) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb 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(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 e(2,3) = e(2,3) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb