Statistiques
| Révision :

root / src / Extrapol_mixed.f90 @ 7

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

1
SUBROUTINE Extrapol_mixed(s,dist,x0,y0,z0,xgeom,Coef)
2

    
3
! This subroutine constructs the path, and if dist<>Infinity, it samples
4
! the path to obtain geometries.
5
! Basically, you call it twice: i) dist=infinity, it will calculate the 
6
! length of the path
7
! ii) dist finite, it will give you the images you want along the path.
8
!
9
! For now, it gives equidistant geometries
10
!!!!!!!!!!!!!!!!
11
!
12
! v2.0
13
! Uses Mix2cart for conversions.
14

    
15

    
16
  use Path_module
17
  use Io_module
18

    
19

    
20
  IMPLICIT NONE
21

    
22

    
23
  REAL(KREAL), INTENT(OUT) :: s
24
  REAL(KREAL), INTENT(IN) :: dist,X0(Nat),Y0(Nat),Z0(Nat)
25
  REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord)
26

    
27
  INTEGER(KINT) :: IdxGeom, I, J, K, Idx
28
  REAL(KREAL) :: Rmsd, MRot(3, 3), ds, u
29

    
30
  REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)
31
  REAL(KREAL), ALLOCATABLE :: IntCoordTmp(:),DerInt(:) ! (NCoord)
32

    
33

    
34
  LOGICAL ::  debug, print
35
  LOGICAL, EXTERNAL :: valid
36

    
37
  !We will calculate the length of the path, in MW coordinates...
38
  ! this is done is a stupid way: we interpolate the zmatrix values,
39
  ! convert them into cartesian, weight the cartesian
40
  ! and calculate the evolution of the distance ! 
41
  ! We have to follow the same procedure for every geometry
42
  ! so even for the first one, we have to convert it from zmat
43
  ! to cartesian !
44

    
45
  debug=valid("pathcreate")
46
  print=valid("printgeom")
47

    
48
  if (debug) Call Header("Entering Extrapol_mixed")
49
  if (debug) WRITE(*,*) "NFroz,NCart",NFroz,NCart
50
  ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),IntCoordTmp(NCoord),DerInt(NCoord))
51
  IdxGeom=1
52
  
53
  IntCoordTmp=IntCoordI(1,1:NCoord)
54
  call Mixed2Cart(Nat,IndZmat,IntCoordTmp,XyzTmp2(1,1))
55

    
56
  XyzTmp=XyzTmp2
57

    
58
! We align this geometry with the original one
59
! PFL 01/feb/2007: useless when dealing with more than three cartesian atoms 
60
! PFL 17/July/2006: only if we have more than 4 atoms.
61
  IF ((NCart.LT.3).AND.(Nat.GT.4)) THEN
62
!  IF (Nat.GT.4) THEN
63
     Call  CalcRmsd(Nat, x0,y0,z0,                              &
64
          xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
65
          MRot,rmsd,FRot,FAlign)
66
  END IF
67

    
68
  XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
69
  IntCoordF(IdxGeom,:)=IntCoordI(1,:)
70

    
71
  if (print.AND.(Dist.LE.1e20)) THEN
72
     WRITE(IOOUT,'(1X,I5)') Nat
73
     WRITE(IOOUT,*) "# Cartesian Coordinates for geom",IdxGeom
74
     DO i=1,Nat
75
        If (Renum) THEN
76
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),      &
77
                (XyzTmp2(Order(I),J),J=1,3)
78
        ELSE
79
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),      &
80
                (XyzTmp2(I,J),J=1,3)
81
        END IF
82
     END DO
83
  END IF
84

    
85
! Calculate tangents for first geometry
86
  u=0.d0
87
     DO Idx=1,NCoord
88
           call splintder(u,IntCoordTmp(Idx),DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
89
     END DO
90
     IntTangent(IdxGeom,:)=DerInt
91

    
92

    
93
! First geometry already initialized
94

    
95
  s=0.
96
!  valzmat=0.d0
97

    
98
  DO K=1,NMaxPtPath
99
     u=real(K)/NMaxPtPath*(NGeomI-1.)
100

    
101
     XYZTmp2=0.
102

    
103
! We generate the interpolated  coordinates
104
     DO Idx=1,NCoord
105
           call splintder(u,IntCoordTmp(Idx),DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
106
     END DO
107

    
108
! We convert it into Cartesian coordinates
109
  call Mixed2Cart(Nat,IndZmat,IntCoordTmp,XyzTmp2(1,1))
110

    
111
! We calculate ds
112
     ds=0.
113
     DO I=1,Nat
114
        DO J=1,3
115
           ds=ds+MassAt(I)*(XYZTMp2(I,J)-XYZTmp(I,J))**2
116
           XYZTmp(I,J)=XyzTMP2(I,J)
117
        ENDDO
118
     ENDDO
119
     s=s+sqrt(ds)
120

    
121
!         if (debug) WRITE(*,*) "Debug u,s,dist",u,s,dist
122
     if (s>=dist) THEN
123
        s=s-dist
124
        IdxGeom=IdxGeom+1
125
        SGeom(IdxGeom)=s+IdxGeom*dist
126
        XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
127

    
128
        IntCoordF(IdxGeom,:)=IntCoordTmp
129
        IntTangent(IdxGeom,:)=DerInt
130

    
131
        if (print) THEN
132
           WRITE(IOOUT,'(1X,A,12(1X,F10.5))') "#Internal Coord",IntCoordTmp(1:NCoord)
133
           WRITE(IOOUT,'(1X,A,12(1X,F10.5))') "#Internal tan",IntTangent(IdxGeom,1:NCoord)
134
           WRITE(IOOUT,'(1X,I5)') Nat
135
           WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K,u,Xgeom(NGeomI)
136

    
137
! PFL 01/feb/2007: useless when dealing with more than three cartesian atoms 
138
! PFL 17/July/2006: only if we have more than 4 atoms.
139
           IF ((NCart.LT.3).AND.(Nat.GT.4)) THEN
140
!           IF (Nat.GT.4) THEN
141
              Call  CalcRmsd(Nat, x0,y0,z0,                               &
142
                   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),   &
143
                   MRot,rmsd,FRot,FAlign)
144
           END IF
145

    
146
           DO i=1,Nat
147
              If (Renum) THEN
148
                 WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),   &
149
                      (XyzTmp2(Order(I),J),J=1,3)
150
              ELSE
151
                 WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),    &
152
                      (XyzTmp2(I,J),J=1,3)
153
              END IF
154
           END DO
155

    
156
        END IF
157
     END IF
158
  END DO
159

    
160

    
161
  if ((s>=0.1*dist).AND.(IdxGeom.EQ.NGeomF)) THEN
162
     WRITE(*,*) "** PathCreate ***"
163
     WRITE(*,*) "Distribution of points along the path is wrong."
164
     WRITE(*,*) "Increase value of NMaxPtPath in the input file"
165
     WRITE(*,*) "Present value is:", NMaxPtPath
166
     STOP
167
  END IF
168

    
169
  IdxGeom=NGeomF
170

    
171
! We have to add the last geometry. We copy the last geom of Initial geometries.
172
  IntCoordF(IdxGeom,:)=IntCoordI(NGeomI,:)
173
  IntTangent(IdxGeom,:)=DerInt
174

    
175
! we convert it to cartesian geom
176
  call Mixed2Cart(Nat,IndZmat,IntCoordF(IdxGeom,1),XyzTmp2(1,1))
177
  XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
178

    
179
     if (print) THEN
180
        WRITE(IOOUT,'(1X,I5)') Nat
181
        WRITE(IOOUT,*) "# Cartesian coord for last Geometry ",IdxGeom,K,u,Xgeom(NGeomI)
182
! PFL 01/feb/2007: useless when dealing with more than three cartesian atoms 
183
! PFL 17/July/2006: only if we have more than 4 atoms.
184
        IF ((NCART.LT.3).AND.(Nat.GT.4)) THEN
185
!       IF (Nat.GT.4) THEN
186
           Call  CalcRmsd(Nat, x0,y0,z0,                          &
187
                xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),    &
188
                MRot,rmsd,FRot,FAlign)
189
        END IF
190

    
191
        DO i=1,Nat
192
           If (Renum) THEN
193
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &
194
                   (XyzTmp2(Order(I),J),J=1,3)
195
           ELSE
196
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),     &
197
                   (XyzTmp2(I,J),J=1,3)
198
           END IF
199
        END DO
200
     END IF
201

    
202

    
203
  if (debug) WRITE(*,*) 's final =',s
204

    
205
  DEALLOCATE(XyzTmp,XyzTmp2)
206

    
207
  if (debug) Call Header("Extrapol_Mixed Over")
208

    
209
 
210
END SUBROUTINE EXTRAPOL_MIXED