Statistiques
| Révision :

root / src / AnaPath.f90 @ 8

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

1
!================================================================
2
!
3
! This subroutine analyzes the geometries along a path
4
! and prints it to FileUnit
5
!================================================================
6

    
7
SUBROUTINE AnaPath(IOpt,FileUnit)
8

    
9

    
10
  use Io_module
11
  use Path_module, only : Nat,  NGeomF, XyzGeomF, SGeom,MassAt, &
12
       IReparam,NbVar, FormAna, Energies,Order, &
13
       Renum, PrintGeomFactor, &
14
       Align, NFroz,  FrozAtoms, AnaFile,NbVar
15

    
16

    
17
  IMPLICIT NONE
18

    
19
  ! Input
20
  ! Iteration number
21
  INTEGER(KINT), INTENT(IN) :: IOpt
22
  ! Unit file to print
23
  INTEGER(KINT), INTENT(IN) :: FileUnit
24

    
25

    
26
  INTERFACE
27
     function valid(string) result (isValid)
28
       CHARACTER(*), intent(in) :: string
29
       logical                  :: isValid
30
     END function VALID
31

    
32
     FUNCTION test_zmat(na,ind_zmat)
33

    
34
       USE Path_module
35

    
36
       logical :: test_zmat
37
       integer(KINT) :: na
38
       integer(KINT) :: ind_zmat(Na,5)
39
     END FUNCTION TEST_ZMAT
40

    
41

    
42
     SUBROUTINE die(routine, msg, file, line, unit)
43

    
44
       Use VarTypes
45
       Use io_module
46

    
47
       implicit none
48

    
49
       character(len=*), intent(in)           :: routine, msg
50
       character(len=*), intent(in), optional :: file
51
       integer(KINT), intent(in), optional      :: line, unit
52

    
53
     END SUBROUTINE die
54

    
55
     SUBROUTINE Warning(routine, msg, file, line, unit)
56

    
57
       Use VarTypes
58
       Use io_module
59

    
60
       implicit none
61

    
62
       character(len=*), intent(in)           :: routine, msg
63
       character(len=*), intent(in), optional :: file
64
       integer(KINT), intent(in), optional      :: line, unit
65

    
66
     END SUBROUTINE Warning
67

    
68
  END INTERFACE
69

    
70
  LOGICAL :: Debug, CalcDist
71
  INTEGER(KINT) :: I,J,K,Iat
72
  REAL(KREAL), ALLOCATABLE :: GeomCart(:,:),GeomCartPrec(:,:) ! Nat,3
73
  REAL(KREAL), ALLOCATABLE :: xRef(:),yRef(:),zRef(:) ! Nat
74
  REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar
75
  REAL(KREAL) ::  ds
76
  REAL(KREAL) :: MRot(3,3), Rmsd  
77
  CHARACTER(SCHARS) :: Plot="plot"
78

    
79
  debug=valid("AnaPath")
80

    
81
  if (debug) Call header("Entering AnaPath")
82

    
83
  !  if (debug) WRITE(*,*) "DBG AnaPaht PrintGeomFactor:",PrintGeomFactor
84
  ALLOCATE(GeomCart(Nat,3))
85
  If (NbVar>0) THEN
86
     ALLOCATE(Values(NbVar))
87
  END IF
88

    
89
  CalcDist=.FALSE.
90

    
91
  WRITE(FileUnit,'("# Iteration ",I5)') Iopt
92

    
93
  DO J=1,Nat
94
     Iat=J
95
     IF (renum) Iat=Order(J)
96
     GeomCart(J,:)=XyzGeomF(1,:,Iat)
97
  END DO
98

    
99

    
100
  ! Do we know the curvilinear values ?
101
  IF (MOD(IOpt,IReparam)/=0) THEN
102
     CalcDist=.TRUE.
103
     SGeom=0.
104
     ALLOCATE(xRef(Nat),yRef(Nat),zref(Nat),GeomCartPrec(Nat,3))
105
     xRef=GeomCart(:,1)
106
     yRef=GeomCart(:,2)
107
     zRef=GeomCart(:,3)
108
     GeomCartPrec=GeomCart
109
  END IF
110

    
111
  if (debug) THEN
112
     WRITE(*,*) "AnaPath  - GeomCart - I=",1
113
     DO K=1,Nat
114
        WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)
115
     END DO
116
     WRITE(*,*) "AnaPath  - XYzGeomF(I,:,:) - I=",1
117
     DO K=1,Nat
118
        WRITE(*,'(1X,I5,3(1X,F15.8))') K,XyzGeomF(1,:,K)
119
     END DO
120
  END IF
121

    
122
  If (NbVar>0) THEN
123
     Call AnalyzeGeom(GeomCart,Values)
124
     WRITE(FileUnit,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1)
125
     if (debug) WRITE(*,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1)
126
  ELSE
127
     WRITE(FileUnit,FormAna) SGeom(1),0.d0,Energies(1)
128
     if (debug) WRITE(*,FormAna) SGeom(1),0.d0,Energies(1)
129
  END IF
130

    
131
  DO I=2,NGeomF
132
     DO J=1,Nat
133
        Iat=J
134
        IF (renum) Iat=Order(J)
135
        GeomCart(J,:)=XyzGeomF(I,:,Iat)
136
     END DO
137

    
138
     If (CalcDist) THEN
139
        ! PFL 2013 Feb 26
140
        ! For now, I do _NOT_ align the geometries
141
        if (debug) THEN
142
           WRITE(*,*) "AnaPath  - GeomCart avant align - I=",I
143
           DO K=1,Nat
144
              WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)
145
           END DO
146
        END IF
147

    
148
        IF ((Nat.GE.4).OR.Align) THEN
149
           if (debug) WRITE(*,*) "DBG AnaPath: Aglin 2 geoms"
150

    
151
           ! If we have frozen atoms we align only those ones.
152
           if (NFroz.GT.0) THEN
153
              Call AlignPartial(Nat,xRef,yRef,zRef,                 &
154
                   GeomCart(1,1),GeomCart(1,2),GeomCart(1,3),       &
155
                   FrozAtoms,MRot)
156
           ELSE 
157
              Call  CalcRmsd(Nat, xRef,yRef,zRef,               &
158
                   GeomCart(1,1),GeomCart(1,2),GeomCart(1,3),   &
159
                   MRot,rmsd,.TRUE.,.TRUE.)
160
           END IF
161
           ! <- PFL 24 Nov 2008
162

    
163
        END IF
164

    
165
        if (debug) WRITE(*,*) "Mass:",MassAt(1:Nat)
166
        ds=0.
167
        DO J=1,Nat
168
           Iat=J
169
           !        IF (renum) Iat=Order(J)
170
           if (debug) WRITE(*,*) "Mass(J):",J,MassAt(Iat)
171

    
172
           do K=1,3
173
              ds=ds+MassAt(Iat)*(GeomCartPrec(J,K)-GeomCart(J,K))**2
174
           END DO
175
           if (debug) WRITE(*,*) "DBG DS:",J,ds
176
        END DO
177
        ds=sqrt(ds) 
178
        if (debug) THEN
179
           WRITE(*,*) "AnaPath  - GeomCartPrec - I=",I
180
           DO K=1,Nat
181
              WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCartPrec(K,1:3)
182
           END DO
183
           WRITE(*,*) "DBG Anapaht, ds=",ds
184
        END IF
185

    
186
        SGeom(I)=SGeom(I-1)+ds
187
        GeomCartPrec=GeomCart
188
        xRef=GeomCart(:,1)
189
        yRef=GeomCart(:,2)
190
        zRef=GeomCart(:,3)
191

    
192
     END IF
193

    
194
     if (debug) THEN
195
        WRITE(*,*) "AnaPath  - GeomCart - I=",I
196
        DO K=1,Nat
197
           WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)
198
        END DO
199
     END IF
200

    
201
     If (NbVar>0) THEN
202
        Call AnalyzeGeom(GeomCart,Values)
203
        WRITE(FileUnit,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*ConvE,Energies(i)
204
        IF (Debug) WRITE(*,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*ConvE,Energies(i)
205
     ELSE
206
        WRITE(FileUnit,FormAna) SGeom(i),(Energies(i)-Energies(1))*ConvE,Energies(i)
207
        If (debug) WRITE(*,FormAna) SGeom(i),(Energies(i)-Energies(1))*ConvE,Energies(i)
208
     END IF
209
  END DO
210

    
211
  WRITE(FileUnit,*) ""
212
  WRITE(FileUnit,*) ""
213

    
214
  DeAllocate(GeomCart)
215
  If (NbVar>0) DeAllocate(Values)
216
  If (CalcDist) Deallocate(GeomCartPrec,xRef,yRef,zRef)
217

    
218
  Write(IoGplot,'(A,I5,A)') TRIM(Plot) // ' "' // TRIM(AnaFile) // '"  i ',Iopt," u xcol:Ecol w lp "
219
  Plot="replot"
220

    
221
  if (debug) Call header("AnaPath over")
222

    
223
END SUBROUTINE AnaPath
224