Statistiques
| Révision :

root / src / CalcTangent.f90 @ 1

Historique | Voir | Annoter | Télécharger (6 ko)

1 1 pfleura2
!This programs reads a set of geometries, and print the tangent
2 1 pfleura2
! to the path they form
3 1 pfleura2
4 1 pfleura2
SUBROUTINE CalcTangent(XgeomF,NGeomS,XGeom,Coef)
5 1 pfleura2
6 1 pfleura2
  use Path_module, only: IntCoordI, Nat, Coord, XyzGeomI, Atome, Order, &
7 1 pfleura2
                         OrderInv, IndZmat, NMaxPtPath, IntTangent, XyzTangent, &
8 1 pfleura2
             NGeomF, NCoord
9 1 pfleura2
  use Io_module, ONly : IoTmp,PathName, Kint,KReal
10 1 pfleura2
  use VarTypes
11 1 pfleura2
12 1 pfleura2
  IMPLICIT NONE
13 1 pfleura2
  ! NGeomS: number of geometries.
14 1 pfleura2
  INTEGER(KINT), INTENT(IN) :: NGeomS
15 1 pfleura2
  REAL(KREAL), INTENT(IN) :: XGeomF(NGeomF),Xgeom(NGeomS),Coef(NGeomS,NCoord)
16 1 pfleura2
17 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: Path(:,:) ! NGeomS,NCoord
18 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: PathTan(:,:) ! NGeomF,NCoord
19 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: DbgTgt(:),GeomTmp(:) ! NCoord
20 1 pfleura2
21 1 pfleura2
  INTEGER(KINT) :: i, j, n3at, Idx, iat,igeom,k
22 1 pfleura2
  INTEGER(KINT), SAVE :: icall=-1
23 1 pfleura2
  INTEGER(KINT) :: NMaxPt
24 1 pfleura2
  REAL(KREAL), PARAMETER :: Inf=1e32, h=0.001d0
25 1 pfleura2
  CHARACTER(120) :: Line,File
26 1 pfleura2
  CHARACTER(LCHARS) :: TITLE
27 1 pfleura2
  LOGICAL :: Debug=.False.
28 1 pfleura2
  REAL(KREAL) :: utmp,der
29 1 pfleura2
30 1 pfleura2
  LOGICAL, EXTERNAL :: valid
31 1 pfleura2
32 1 pfleura2
  debug=valid('calctangent')
33 1 pfleura2
34 1 pfleura2
  if (debug) WRITE(*,*) "======================================== Entering CalcTangen ============================"
35 1 pfleura2
  n3at=3*nat
36 1 pfleura2
37 1 pfleura2
  ALLOCATE(Path(NGeomS,Ncoord),PathTan(NGeomF,NCoord),GeomTMP(NCoord))
38 1 pfleura2
39 1 pfleura2
 icall=icall+2
40 1 pfleura2
 write(Line,'(I5)') icall
41 1 pfleura2
 File=TRIM(adjustl(PathName)) // '_dbgtgt_' // TRIM(adjustl(Line)) // '.dat'
42 1 pfleura2
 OPEN(IOTMP,File=File)
43 1 pfleura2
44 1 pfleura2
  ! We construct path, depending on COORD:
45 1 pfleura2
  ! Path(NGeomS,Ncoord), NGeoms is equal to NGeomI
46 1 pfleura2
  ! IntCoordI(:,:) ! (NGeomI,3*Nat-6)
47 1 pfleura2
 SELECT CASE(Coord)
48 1 pfleura2
    CASE ('ZMAT','MIXED')
49 1 pfleura2
     Path(:,:)=IntCoordI(:,:)
50 1 pfleura2
  CASE ('CART','HYBRID')
51 1 pfleura2
     Path=reshape(XyzGeomI,(/NGeomS,NCoord/))
52 1 pfleura2
  CASE ('BAKER')
53 1 pfleura2
   Path(:,:)=IntCoordI(:,:)
54 1 pfleura2
  CASE DEFAULT
55 1 pfleura2
     WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in CalcTangent. STOP"
56 1 pfleura2
     STOP
57 1 pfleura2
  END SELECT
58 1 pfleura2
59 1 pfleura2
60 1 pfleura2
  if (debug) THEN
61 1 pfleura2
     WRITE(*,*) 'DBG CalcTangent. Geometries from Path'
62 1 pfleura2
     DO I=1,NGeomS
63 1 pfleura2
        GeomTmp=Path(I,1:NCoord)
64 1 pfleura2
        WRITE(Title,"('Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomS
65 1 pfleura2
        Call PrintGeom(Title,Nat,Ncoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
66 1 pfleura2
!         IF (COORD.EQ.'ZMAT') THEN
67 1 pfleura2
!            WRITE(*,"('Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomS
68 1 pfleura2
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
69 1 pfleura2
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
70 1 pfleura2
!                 Nom(Atome(OrderInv(indzmat(2,2)))),Path(I,1)
71 1 pfleura2
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
72 1 pfleura2
!                 Nom(Atome(OrderInv(indzmat(3,2)))),Path(I,2), &
73 1 pfleura2
!                 Nom(Atome(OrderInv(indzmat(3,3)))),Path(I,3)
74 1 pfleura2
!            Idx=4
75 1 pfleura2
!            DO Iat=4,Nat
76 1 pfleura2
!               WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
77 1 pfleura2
!                    Nom(Atome(OrderInv(indzmat(iat,2)))),Path(I,Idx), &
78 1 pfleura2
!                    Nom(Atome(OrderInv(indzmat(iat,3)))),Path(I,Idx+1), &
79 1 pfleura2
!                    Nom(Atome(OrderInv(indzmat(iat,4)))),Path(I,Idx+2)
80 1 pfleura2
!               Idx=Idx+3
81 1 pfleura2
!            END DO
82 1 pfleura2
!         ELSE
83 1 pfleura2
!            WRITE(*,*) Nat
84 1 pfleura2
!            WRITE(*,"('Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF
85 1 pfleura2
!            DO Iat=1,Nat
86 1 pfleura2
!               WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
87 1 pfleura2
!                    Path(I,3*Iat-2:3*Iat)
88 1 pfleura2
!            END DO
89 1 pfleura2
!         END IF
90 1 pfleura2
     END DO
91 1 pfleura2
  END IF
92 1 pfleura2
93 1 pfleura2
!Debug info
94 1 pfleura2
 IF (debug.AND.(COORD.EQ.'ZMAT')) THEN
95 1 pfleura2
     ALLOCATE(DbgTgt(NCoord))
96 1 pfleura2
     WRITE(*,*) 'DBG CalcTangent. Path written in ' // TRIM(File)
97 1 pfleura2
     WRITE(*,'(A,12(1X,F10.5))') "DBG CalcTangent, XgeomF:",XGeomF
98 1 pfleura2
99 1 pfleura2
     NMaxPt=min(200,NMaxPtPath)-1
100 1 pfleura2
    DO K=0,NMaxPt
101 1 pfleura2
       utmp=K*Xgeom(NGeomS)/NMaxPt
102 1 pfleura2
       DO J=1,NCoord
103 1 pfleura2
       Call Splint(utmp,DbgTgt(J),NGeomS,Xgeom,Path(1,J),Coef(1,J))
104 1 pfleura2
       END DO
105 1 pfleura2
       WRITE(IOTMP,'(12(1X,F10.5))') utmp,DbgTgt
106 1 pfleura2
    END DO
107 1 pfleura2
    DEALLOCATE(DbgTgt)
108 1 pfleura2
 END IF
109 1 pfleura2
110 1 pfleura2
!We construct the derivatives
111 1 pfleura2
DO I=1,NGeomF
112 1 pfleura2
  DO J=1,NCoord
113 1 pfleura2
        Call SplinDer(XgeomF(I),PathTan(I,J),NGeomS,Xgeom,Path(1,J),Coef(1,J))
114 1 pfleura2
  END DO
115 1 pfleura2
END DO
116 1 pfleura2
117 1 pfleura2
 SELECT CASE(Coord)
118 1 pfleura2
    CASE ('ZMAT','MIXED')
119 1 pfleura2
     IntTangent=PathTan
120 1 pfleura2
  CASE ('CART','HYBRID')
121 1 pfleura2
     XyzTangent=Pathtan
122 1 pfleura2
  END SELECT
123 1 pfleura2
124 1 pfleura2
125 1 pfleura2
  if (debug) THEN
126 1 pfleura2
     WRITE(*,*) 'DBG CalcTangent. Tangents'
127 1 pfleura2
     DO I=1,NGeomF
128 1 pfleura2
        GeomTmp=PathTan(I,1:NCoord)
129 1 pfleura2
        WRITE(Title,"('Tangents for geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF
130 1 pfleura2
        Call PrintGeom(Title,Nat,Ncoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
131 1 pfleura2
132 1 pfleura2
        WRITe(*,'(A,12(1X,F10.6))') "Test:",GeomTmp(1:NCoord)
133 1 pfleura2
134 1 pfleura2
!         IF (COORD.EQ.'ZMAT') THEN
135 1 pfleura2
!            WRITE(*,"('Tangent for Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF
136 1 pfleura2
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
137 1 pfleura2
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
138 1 pfleura2
!                 Nom(Atome(OrderInv(indzmat(2,2)))),IntTangent(I,1)
139 1 pfleura2
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
140 1 pfleura2
!                 Nom(Atome(OrderInv(indzmat(3,2)))),IntTangent(I,2), &
141 1 pfleura2
!                 Nom(Atome(OrderInv(indzmat(3,3)))),IntTangent(I,3)
142 1 pfleura2
!            Idx=4
143 1 pfleura2
!            DO Iat=4,Nat
144 1 pfleura2
!               WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
145 1 pfleura2
!                    Nom(Atome(OrderInv(indzmat(iat,2)))),IntTangent(I,Idx), &
146 1 pfleura2
!                    Nom(Atome(OrderInv(indzmat(iat,3)))),IntTangent(I,Idx+1), &
147 1 pfleura2
!                    Nom(Atome(OrderInv(indzmat(iat,4)))),IntTangent(I,Idx+2)
148 1 pfleura2
!               Idx=Idx+3
149 1 pfleura2
!            END DO
150 1 pfleura2
!         ELSE
151 1 pfleura2
!            WRITE(*,*) Nat
152 1 pfleura2
!            WRITE(*,"('Tangent for Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF
153 1 pfleura2
!            DO Iat=1,Nat
154 1 pfleura2
!               WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
155 1 pfleura2
!                    XyzTangent(I,3*Iat-2:3*Iat)
156 1 pfleura2
!            END DO
157 1 pfleura2
!         END IF
158 1 pfleura2
     END DO
159 1 pfleura2
  END IF
160 1 pfleura2
161 1 pfleura2
162 1 pfleura2
DeALLOCATE(Path)
163 1 pfleura2
164 1 pfleura2
 CLOSE(IOTMP)
165 1 pfleura2
  if (debug) WRITE(*,*) "======================================== Ending CalcTangent ============================"
166 1 pfleura2
167 1 pfleura2
END SUBROUTINE CalcTangent