Statistiques
| Révision :

root / src / AnaPath.f90 @ 7

Historique | Voir | Annoter | Télécharger (5,31 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, AtName,SGeom,MassAt, &
12
       IReparam,NbVar, FormAna, Energies, au2kcal,XyzGeomI,Order,OrderInv, &
13
       Renum, PrintGeomFactor, &
14
       Align, NFroz,  FrozAtoms
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

    
78
  debug=valid("AnaPath")
79

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

    
82
  if (debug) WRITE(*,*) "DBG AnaPaht PrintGeomFactor:",PrintGeomFactor
83
  ALLOCATE(GeomCart(Nat,3))
84
  ALLOCATE(Values(NbVar))
85

    
86
  CalcDist=.FALSE.
87

    
88
  WRITE(FileUnit,'("# Iteration ",I5)') Iopt
89

    
90
  DO J=1,Nat
91
     Iat=J
92
     IF (renum) Iat=Order(J)
93
     GeomCart(J,:)=XyzGeomF(1,:,Iat)
94
  END DO
95

    
96

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

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

    
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

    
123
  DO I=2,NGeomF
124
     DO J=1,Nat
125
        Iat=J
126
        IF (renum) Iat=Order(J)
127
        GeomCart(J,:)=XyzGeomF(I,:,Iat)
128
     END DO
129

    
130
     If (CalcDist) THEN
131
! PFL 2013 Feb 26
132
! For now, I do _NOT_ align the geometries
133
        if (debug) THEN
134
           WRITE(*,*) "AnaPath  - GeomCart avant align - I=",I
135
           DO K=1,Nat
136
              WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)
137
           END DO
138
        END IF
139

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

    
143
! If we have frozen atoms we align only those ones.
144
           if (NFroz.GT.0) THEN
145
              Call AlignPartial(Nat,xRef,yRef,zRef,                 &
146
                   GeomCart(1,1),GeomCart(1,2),GeomCart(1,3),       &
147
                   FrozAtoms,MRot)
148
           ELSE 
149
              Call  CalcRmsd(Nat, xRef,yRef,zRef,               &
150
                   GeomCart(1,1),GeomCart(1,2),GeomCart(1,3),   &
151
                   MRot,rmsd,.TRUE.,.TRUE.)
152
           END IF
153
              ! <- PFL 24 Nov 2008
154
           
155
        END IF
156

    
157
        if (debug) WRITE(*,*) "Mass:",MassAt(1:Nat)
158
        ds=0.
159
        DO J=1,Nat
160
           Iat=J
161
!        IF (renum) Iat=Order(J)
162
           if (debug) WRITE(*,*) "Mass(J):",J,MassAt(Iat)
163
            
164
           do K=1,3
165
              ds=ds+MassAt(Iat)*(GeomCartPrec(J,K)-GeomCart(J,K))**2
166
           END DO
167
           if (debug) WRITE(*,*) "DBG DS:",J,ds
168
        END DO
169
        ds=sqrt(ds)
170
        if (debug) THEN
171
           WRITE(*,*) "AnaPath  - GeomCartPrec - I=",I
172
           DO K=1,Nat
173
              WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCartPrec(K,1:3)
174
           END DO
175
           WRITE(*,*) "DBG Anapaht, ds=",ds
176
        END IF
177

    
178
        SGeom(I)=SGeom(I-1)+ds
179
        GeomCartPrec=GeomCart
180
        xRef=GeomCart(:,1)
181
        yRef=GeomCart(:,2)
182
        zRef=GeomCart(:,3)
183
        
184
     END IF
185

    
186
     if (debug) THEN
187
        WRITE(*,*) "AnaPath  - GeomCart - I=",I
188
        DO K=1,Nat
189
           WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)
190
        END DO
191
     END IF
192

    
193
     Call AnalyzeGeom(GeomCart,Values)
194
     WRITE(FileUnit,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*au2kcal,Energies(i)
195
  END DO
196

    
197
  WRITE(FileUnit,*) ""
198
  WRITE(FileUnit,*) ""
199

    
200
  DeAllocate(GeomCart,Values)
201
  If (CalcDist) Deallocate(GeomCartPrec,xRef,yRef,zRef)
202

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

    
205
END SUBROUTINE AnaPath
206