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 |
|