fix the model and nnadia

This commit is contained in:
jean paul nshuti 2025-10-09 13:01:49 +02:00
parent ccf300a77e
commit 6f1323639a
7 changed files with 1295 additions and 494 deletions

View File

@ -0,0 +1,869 @@
Module dip_param
IMPLICIT NONE
Integer,parameter :: np=58
Double precision :: p(58)
integer :: pst(2,400)
contains
SUBROUTINE init_dip_planar_data()
implicit none
p( 1)= 0.101430388D+01
p( 2)= 0.000000000D+00
p( 3)= 0.157490804D+01
p( 4)= 0.000000000D+00
p( 5)= -0.343280555D+01
p( 6)= 0.000000000D+00
p( 7)= 0.246060523D+00
p( 8)= 0.000000000D+00
p( 9)= 0.000000000D+00
p( 10)= 0.408099282D-01
p( 11)= 0.000000000D+00
p( 12)= 0.000000000D+00
p( 13)= -0.769396889D+01
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.000000000D+00
p( 22)= 0.217421464D+00
p( 23)= 0.464794515D+00
p( 24)= 0.000000000D+00
p( 25)= -0.243354772D+00
p( 26)= 0.000000000D+00
p( 27)= 0.000000000D+00
p( 28)= 0.000000000D+00
p( 29)= 0.000000000D+00
p( 30)= 0.000000000D+00
p( 31)= 0.000000000D+00
p( 32)= 0.000000000D+00
p( 33)= 0.000000000D+00
p( 34)= 0.000000000D+00
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.000000000D+00
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.000000000D+00
p( 50)= 0.000000000D+00
p( 51)= 0.000000000D+00
p( 52)= -0.565000000D-02
p( 53)= 0.180418089D+00
p( 54)= 0.000000000D+00
p( 55)= 0.118050680D-01
p( 56)= 0.000000000D+00
p( 57)= 0.000000000D+00
p( 58)= 0.000000000D+00
pst(1, 1)= 1
pst(2, 1)= 2
pst(1, 2)= 3
pst(2, 2)= 2
pst(1, 3)= 5
pst(2, 3)= 2
pst(1, 4)= 7
pst(2, 4)= 3
pst(1, 5)= 10
pst(2, 5)= 3
pst(1, 6)= 13
pst(2, 6)= 3
pst(1, 7)= 16
pst(2, 7)= 2
pst(1, 8)= 18
pst(2, 8)= 2
pst(1, 9)= 20
pst(2, 9)= 2
pst(1, 10)= 22
pst(2, 10)= 1
pst(1, 11)= 23
pst(2, 11)= 2
pst(1, 12)= 25
pst(2, 12)= 5
pst(1, 13)= 30
pst(2, 13)= 10
pst(1, 14)= 40
pst(2, 14)= 1
pst(1, 15)= 41
pst(2, 15)= 2
pst(1, 16)= 43
pst(2, 16)= 4
pst(1, 17)= 47
pst(2, 17)= 0
pst(1, 18)= 47
pst(2, 18)= 2
pst(1, 19)= 49
pst(2, 19)= 3
pst(1, 20)= 52
pst(2, 20)= 1
pst(1, 21)= 53
pst(2, 21)= 2
pst(1, 22)= 55
pst(2, 22)= 4
pst(1, 23)= 0
pst(2, 23)= 0
pst(1, 24)= 0
pst(2, 24)= 0
pst(1, 25)= 0
pst(2, 25)= 0
pst(1, 26)= 0
pst(2, 26)= 0
pst(1, 27)= 0
pst(2, 27)= 0
pst(1, 28)= 0
pst(2, 28)= 0
pst(1, 29)= 0
pst(2, 29)= 0
pst(1, 30)= 0
pst(2, 30)= 0
pst(1, 31)= 0
pst(2, 31)= 0
pst(1, 32)= 0
pst(2, 32)= 0
pst(1, 33)= 0
pst(2, 33)= 0
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

View File

@ -1 +0,0 @@
../../../../Genetic/NO3/Dipole_NO3/Fit_stretch_Latest/fit_genric_bend_no3.f90

396
src/model/model_nh3.f90 Normal file
View File

@ -0,0 +1,396 @@
module diabmodel
use iso_fortran_env, only: idp => int32, dp => real64
use dip_param
implicit none
include "nnparams.incl"
integer(idp),parameter:: ndiab=4
logical :: debug=.false.
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(ndiab,ndiab)
integer(idp) id,i,j ,ii
real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
real(dp),dimension(12):: shift,scal
xs=q(5)
ys=q(6)
xb=q(7)
yb=q(8)
a=q(4)
b=q(9)
call init_dip_planar_data()
! modify the parametr
shift=1.0_dp
scal=1.0d-3
! V term of A2''
ii=1
do i =1,np
if (p(i) .ne. 0) then
p(i) =p(i)*(shift(ii) + scal(ii)*nn_out(ii) )
ii=ii+1
else
p(i) = p(i)
endif
enddo
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=1 !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(ndiab,ndiab)
integer(idp) id,i,j, ii
real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
real(dp),dimension(12):: shift,scal
xs=q(5)
ys=q(6)
xb=q(7)
yb=q(8)
a=q(4)
b=q(9)
call init_dip_planar_data()
! modify the parametr
shift=1.0_dp
scal=1.0d-3
! V term of A2''
ii=1
do i =1,np
if (p(i) .ne. 0) then
p(i) =p(i)*(shift(ii) + scal(ii)*nn_out(ii) )
ii=ii+1
else
p(i) = p(i)
endif
enddo
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=1 !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)
end subroutine diab_y
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

View File

@ -1,469 +0,0 @@
module diabmodel
!use dim_parameter,only:qn,ndiab,pst
use iso_fortran_env, only: idp => int32, dp => real64
!use accuracy_constants, only:dp,idp
use dip_param, only: init_dip_planar_data, p,pst,np
implicit none
include "nnparams.incl"
integer(idp),parameter:: ndiab=5
logical :: debug=.false.
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine diab_x(q,e,nn_out)
real(dp),intent(in)::q(maxnin)
real(dp),intent(out)::e(ndiab,ndiab)
real(dp),intent(inout):: nn_out(maxnout)
integer(idp) id,i,j, ii, non_zer_p
real(dp) xs,xb,ys,yb,a,b,ss,sb,v3(8),v2(6)
xs=q(5)
ys=q(6)
xb=q(7)
yb=q(8)
a=q(4)
b=q(9)
call init_dip_planar_data()
non_zer_p = count(p /= 0.0d0)
ii=1
do i=1,np-2
if (p(i) .ne. 0) then
p(i) =p(i)*(1.0_dp + 1.0d-2 + nn_out(ii))
ii=ii+1
else
p(i)=p(i)
endif
enddo
ss=xs**2+ys**2 ! totaly symmetric term
sb=xb**2+yb**2
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.0d0
id=1 ! 1
e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 param
id=id+1 ! 2
e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 p
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 ! 2 p
e(5,5)=e(5,5)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
id=id+1 ! 4
! order 2
e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & ! 5 p
+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)
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)
! 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
! 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))
!e(2,3)=e(2,3)
! order 1
id=id+1 ! 11 ! 2 param
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 ! 3p
do i=1,3
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(2,3)=e(2,3)+ p(pst(1,id)+(i-1))*v2(i+3)
enddo
! 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
!e(2,3)=e(2,3)- p(pst(1,id))*ys*ss +p(pst(1,id)+1)*ss*2*xs*ys
! W and Z for E2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
id=id+1 ! 14
e(4,4)=e(4,4)+p(pst(1,id))
e(5,5)=e(5,5)-p(pst(1,id))
e(4,5)=e(4,5)
! order 1
id=id+1 ! 2 param 15
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(4,5)=e(4,5)- p(pst(1,id))*ys-p(pst(1,id)+1)*yb
! order 2
id=id+1 ! 16 ! 3p
do i=1,3
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(4,5)=e(4,5)+ p(pst(1,id)+(i-1))*v2(i+3)
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
!e(4,4) = e(4,4)+ p(pst(1,id))*xs*ss + p(pst(1,id)+1)*xb*sb
!e(5,5)=e(5,5)- (p(pst(1,id))*xs*ss + p(pst(1,id)+1)*xb*sb)
!e(4,5)=e(4,5)- p(pst(1,id))*ys*ss -p(pst(1,id)+1)*yb*sb
! E1 X E2
! WW and ZZ
id =id+1 ! 18
e(2,4)=e(2,4)+p(pst(1,id))*b
e(3,5)=e(3,5)-p(pst(1,id))*b
! ORDER 1
id=id+1 ! 19 ! 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(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(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
id=id+1 ! 20
do i=1,3 ! 9 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(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(3,4)=e(3,4)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i+3)
enddo
! pseudo A2 & E1
! ##################################################
!###################################################
! order 0
id=id+1 ! 1 param ! 21
e(1,3)=e(1,3)+b*(p(pst(1,id)))
! order 1
id = id +1 ! 22
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 ! 23
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))
! COUPLING OF A2 WITH E2
!##########################################################################################################
! order 0
id =id+1 !24
e(1,5)=e(1,5)+p(pst(1,id))
! order 1
id = id +1 ! 25
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)
! order 2
id=id+1 ! 26
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) &
+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
call copy_2_lower_triangle(e)
!temp = e(2,2)
!e(2,2)=e(3,3)
!e(3,3)=temp
if (debug) then
do i=1,ndiab
write(34,'(5f14.6)') (e(i,j),j=1,ndiab)
enddo
write(34,*)""
endif
end subroutine diab_x
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! THE Y COMPONENT OF DIPOLE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine diab_y(q,e,nn_out)
real(dp),intent(in)::q(maxnin)
real(dp),intent(out)::e(ndiab,ndiab)
real(dp),intent(inout):: nn_out(maxnout)
integer(idp) id,i,j, ii, non_zer_p
real(dp) xs,xb,ys,yb,a,b,ss,sb,v3(8),v2(6)
xs=q(5)
ys=q(6)
xb=q(7)
yb=q(8)
a=q(4)
b=q(9)
call init_dip_planar_data()
non_zer_p = count(p /= 0.0d0)
ii=1
do i=1,np-2
if (p(i) .ne. 0) then
p(i) =p(i)*(1.0_dp + 1.0d-2 + nn_out(ii))
ii=ii+1
else
p(i)=p(i)
endif
enddo
ss=xs**2+ys**2 ! totaly symmetric term
sb=xb**2+yb**2
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.0d0
! V-term
id=1 ! 1
! order 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
e(5,5)=e(5,5)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb
id=id+1 ! 4b*(
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))
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
! order 0
id=id+1 ! 10
e(2,3)=e(2,3)+p(pst(1,id))
! order 1
id=id+1 !
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
! order 2
id=id+1 ! 12
do i=1,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(2,3)=e(2,3)-p(pst(1,id)+(i-1))*v2(i)
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
!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! order 0
id=id+1 ! 14
e(4,5)=e(4,5)+p(pst(1,id))
! order 1
id=id+1 ! 15
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(4,5)=e(4,5)-p(pst(1,id))*xs -p(pst(1,id)+1)*xb
! order 2
id=id+1 ! 16
do i=1,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(4,5)=e(4,5)-p(pst(1,id)+(i-1))*v2(i)
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
!ORDER 0
id=id+1 ! 18
e(2,5)=e(2,5)+p(pst(1,id))*b
e(3,4)=e(3,4)+p(pst(1,id))*b
! order 1
id=id+1
e(2,4)=e(2,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,5)=e(3,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))*xs+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb)
e(3,4)=e(3,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)
! order 2
id=id+1
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)*b
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)*b
e(2,5)=e(2,5)+(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i)*b
e(3,4)=e(3,4)+(-p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i)*b
! no order 3
!!!!!!!!!!!!!!!!
! the coupling A2 & E1
! #####################
! order 0
id=id+1
e(1,2)=e(1,2)+b*(p(pst(1,id)))
! order 1
id=id+1
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
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))
! COUPLING OF A2 WITH E2
!#######################################################################################
!###############################################################################
! order 0
id = id+1
e(1,4)=e(1,4)+p(pst(1,id))
! order 1
id=id+1
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)
! order 2
id=id+1
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))
e(1,5)=e(1,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))
!write(*,*)'idy=',id
call copy_2_lower_triangle(e)
if (debug) then
do i=1,ndiab
write(*,'(5f14.6)') (e(i,j),j=1,ndiab)
enddo
write(*,*)""
endif
end subroutine diab_y
subroutine copy_2_lower_triangle(mat)
real(dp), intent(inout) :: mat(:, :)
integer :: m, n
! write lower triangle of matrix symmetrical
do n=1,size(mat,1)
do m=n,size(mat,1)
mat(m,n)=mat(n,m)
enddo
enddo
end subroutine copy_2_lower_triangle
end module diabmodel

View File

@ -6,7 +6,7 @@
subroutine nnadia(coords,coeffs,adiaoutp) subroutine nnadia(coords,coeffs,adiaoutp)
USE diabmodel, only: diab_x USE diabmodel, only: diab_x, diab_y
implicit none implicit none
! returns THE DIABATIC MATRIX UPPER MATRIX OF BOTH X AND Y COMPONENT OF DIPOLE ! returns THE DIABATIC MATRIX UPPER MATRIX OF BOTH X AND Y COMPONENT OF DIPOLE
! !
@ -23,19 +23,19 @@
! eigenvs: dummy storage for eigenvectors ! eigenvs: dummy storage for eigenvectors
include 'nnparams.incl' include 'nnparams.incl'
include 'JTmod.incl' !include 'JTmod.incl'
integer ndiab integer ndiab
parameter ndiab=5 parameter ndiab=4
DOUBLE PRECISION coords(maxnin),coeffs(maxnout) DOUBLE PRECISION coords(maxnin),coeffs(maxnout)
DOUBLE PRECISION adiaoutp(maxpout),dmat_x(ndiab,ndiab) DOUBLE PRECISION adiaoutp(maxpout),dmat_x(ndiab,ndiab)
!DOUBLE PRECISION dmat_y(ndiab,ndiab) DOUBLE PRECISION dmat_y(ndiab,ndiab)
!integer i,j,ii !integer i,j,ii
call diab_x(coords,dmat_x,coeffs) call diab_x(dmat_x,coords,coeffs)
!call diab_y(coords,dmat_y,coeffs) call diab_y(dmat_y,coords,coeffs)
! rewrite the diabatic matrix into 1D array of adiaoutp ! rewrite the diabatic matrix into 1D array of adiaoutp
!ii=1 !ii=1
!do i=1,ndiab !do i=1,ndiab
!do j=i,ndiab !do j=i,ndiab
! adiaoutp(ii)=dmat_x(i,j) ! adiaoutp(ii)=dmat_x(i,j)
! adiaoutp(ii+10) = -1*dmat_y(i,j) ! adiaoutp(ii+10) = -1*dmat_y(i,j)
@ -43,20 +43,25 @@
!enddo !enddo
!enddo !enddo
adiaoutp(1) = dmat_x(1,1) adiaoutp(1) = dmat_x(1,1)
adiaoutp(2) = dmat_x(2,1) adiaoutp(2) = dmat_x(1,2)
adiaoutp(3) = dmat_x(2,2) adiaoutp(3) = dmat_x(1,3)
adiaoutp(4) = dmat_x(3,1) adiaoutp(4) = dmat_x(1,4)
adiaoutp(5) = dmat_x(3,2) adiaoutp(5) = dmat_x(2,3)
adiaoutp(6) = dmat_x(3,3) adiaoutp(6) = dmat_x(2,2)
adiaoutp(7) = dmat_x(4,1) adiaoutp(7) = dmat_x(3,1)
adiaoutp(8) = dmat_x(4,2) adiaoutp(8) = dmat_x(3,2)
adiaoutp(9) = dmat_x(4,3) adiaoutp(9) = dmat_x(3,3)
adiaoutp(10) = dmat_x(4,4) adiaoutp(10) = dmat_x(4,4)
adiaoutp(11) = dmat_x(5,1) adiaoutp(11) = -1.0*dmat_y(1,1)
adiaoutp(12) = dmat_x(5,2) adiaoutp(12) = -1.0*dmat_y(1,2)
adiaoutp(13) = dmat_x(5,3) adiaoutp(13) = -1.0*dmat_y(1,3)
adiaoutp(14) = dmat_x(5,4) adiaoutp(14) = -1.0*dmat_y(1,4)
adiaoutp(15) = dmat_x(5,5) adiaoutp(15) = -1.0*dmat_y(2,3)
adiaoutp(16) = -1.0*dmat_y(2,2)
adiaoutp(17) = -1.0*dmat_y(3,1)
adiaoutp(18) = -1.0*dmat_y(3,2)
adiaoutp(19) = -1.0*dmat_y(3,3)
adiaoutp(20) = -1.0*dmat_y(4,4)
!write(*,*) dmat_x(1,1) !write(*,*) dmat_x(1,1)
END SUBROUTINE END SUBROUTINE

View File

@ -51,7 +51,7 @@
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
module ctrans_mod module ctrans_mod
implicit none implicit none
include 'only_model.incl' !include 'only_model.incl'
include 'nnparams.incl' include 'nnparams.incl'
! precalculate pi, 2*pi and angle to radian conversion ! precalculate pi, 2*pi and angle to radian conversion
double precision, parameter :: pii = 4.00d0*datan(1.00d0) double precision, parameter :: pii = 4.00d0*datan(1.00d0)
@ -62,7 +62,7 @@ module ctrans_mod
double precision, parameter:: sq3 = 1.00d0/dsqrt(3.00d0) double precision, parameter:: sq3 = 1.00d0/dsqrt(3.00d0)
double precision, parameter:: sq6 = 1.00d0/dsqrt(6.00d0) double precision, parameter:: sq6 = 1.00d0/dsqrt(6.00d0)
! change distances for equilibrium ! change distances for equilibrium
double precision, parameter :: dchequi = 2.344419d0 double precision, parameter :: dchequi = 1.0228710942d0 ! req NH3+
contains contains
subroutine ctrans(q) subroutine ctrans(q)

1
src/model/nnparams.incl Symbolic link
View File

@ -0,0 +1 @@
../nnparams.incl