Révision 8 src/AnaPath.f90

AnaPath.f90 (revision 8)
8 8

  
9 9

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

  
16 16

  
17 17
  IMPLICIT NONE
18 18

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

  
25 25

  
......
39 39
     END FUNCTION TEST_ZMAT
40 40

  
41 41

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

  
44
      Use VarTypes
45
      Use io_module
44
       Use VarTypes
45
       Use io_module
46 46

  
47
      implicit none
47
       implicit none
48 48

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

  
53
    END SUBROUTINE die
53
     END SUBROUTINE die
54 54

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

  
57
      Use VarTypes
58
      Use io_module
57
       Use VarTypes
58
       Use io_module
59 59

  
60
      implicit none
60
       implicit none
61 61

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

  
66
    END SUBROUTINE Warning
66
     END SUBROUTINE Warning
67 67

  
68 68
  END INTERFACE
69 69

  
......
74 74
  REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar
75 75
  REAL(KREAL) ::  ds
76 76
  REAL(KREAL) :: MRot(3,3), Rmsd  
77
  CHARACTER(SCHARS) :: Plot="plot"
77 78

  
78 79
  debug=valid("AnaPath")
79 80

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

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

  
86 89
  CalcDist=.FALSE.
87 90

  
......
94 97
  END DO
95 98

  
96 99

  
97
! Do we know the curvilinear values ?
100
  ! Do we know the curvilinear values ?
98 101
  IF (MOD(IOpt,IReparam)/=0) THEN
99 102
     CalcDist=.TRUE.
100 103
     SGeom=0.
......
116 119
     END DO
117 120
  END IF
118 121

  
119
  Call AnalyzeGeom(GeomCart,Values)
120
  WRITE(FileUnit,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1)
121
  if (debug) WRITE(*,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1)
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
122 130

  
123 131
  DO I=2,NGeomF
124 132
     DO J=1,Nat
......
128 136
     END DO
129 137

  
130 138
     If (CalcDist) THEN
131
! PFL 2013 Feb 26
132
! For now, I do _NOT_ align the geometries
139
        ! PFL 2013 Feb 26
140
        ! For now, I do _NOT_ align the geometries
133 141
        if (debug) THEN
134 142
           WRITE(*,*) "AnaPath  - GeomCart avant align - I=",I
135 143
           DO K=1,Nat
......
138 146
        END IF
139 147

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

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

  
155 163
        END IF
156 164

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

  
164 172
           do K=1,3
165 173
              ds=ds+MassAt(Iat)*(GeomCartPrec(J,K)-GeomCart(J,K))**2
166 174
           END DO
167 175
           if (debug) WRITE(*,*) "DBG DS:",J,ds
168 176
        END DO
169
        ds=sqrt(ds)
177
        ds=sqrt(ds) 
170 178
        if (debug) THEN
171 179
           WRITE(*,*) "AnaPath  - GeomCartPrec - I=",I
172 180
           DO K=1,Nat
......
180 188
        xRef=GeomCart(:,1)
181 189
        yRef=GeomCart(:,2)
182 190
        zRef=GeomCart(:,3)
183
        
191

  
184 192
     END IF
185 193

  
186 194
     if (debug) THEN
......
190 198
        END DO
191 199
     END IF
192 200

  
193
     Call AnalyzeGeom(GeomCart,Values)
194
     WRITE(FileUnit,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*au2kcal,Energies(i)
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
195 209
  END DO
196 210

  
197 211
  WRITE(FileUnit,*) ""
198 212
  WRITE(FileUnit,*) ""
199 213

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

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

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

  
205 223
END SUBROUTINE AnaPath

Formats disponibles : Unified diff