root / src / AnaPath.f90 @ 8
History | View | Annotate | Download (6 kB)
| 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 |
|