root / src / Extrapol_baker.f90 @ 2
Historique | Voir | Annoter | Télécharger (10,67 ko)
| 1 |
SUBROUTINE Extrapol_baker(s,dist,x0,y0,z0,xgeom,Coef,XgeomF) |
|---|---|
| 2 |
|
| 3 |
! This subroutine constructs the path, andabscissa if dist<>Infinity, it samples |
| 4 |
! the path to obtain geometries. |
| 5 |
! Basically, you call it twice: i) dist=infinity, it will calculate the length of the path |
| 6 |
! ii) dist finite, it will give you the images you want along the path. |
| 7 |
! |
| 8 |
! For now, it gives equidistant geometries. |
| 9 |
! |
| 10 |
! A reference geometry for the alignment: X0(Nat),Y0(Nat),Z0(Nat) |
| 11 |
|
| 12 |
use Path_module, only : IntCoordI, NMaxPtPath, XyzGeomF, IntCoordF, & |
| 13 |
IntTangent, Renum, Nom, Order, MassAt, SGeom, Nat, NGeomI, & |
| 14 |
NGeomF, Atome, NCoord, OrderInv, XyzGeomI,BTransInvF, & |
| 15 |
XPrimitive,XPrimitiveF, NPrim, & |
| 16 |
BTransInv_local,UMatF,UMat_local,FirstTimePathCreate,Pi |
| 17 |
! IntCoordI(NGeomI,3*Nat-6), Coef(NGeomI,NCoord), NMaxPtPath=1000 (default value) |
| 18 |
! More appropriate: IntCoordI(NGeomI,NCoord) |
| 19 |
use Io_module |
| 20 |
IMPLICIT NONE |
| 21 |
|
| 22 |
REAL(KREAL), INTENT(OUT) :: s |
| 23 |
! A reference geometry for the alignment: |
| 24 |
REAL(KREAL), INTENT(IN) :: dist,X0(Nat),Y0(Nat),Z0(Nat) |
| 25 |
! Xgeom(NGeomI): abscissa of all initial geometries. |
| 26 |
! Coef(NGeomI,NCoord): spline coefficients. |
| 27 |
REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord) |
| 28 |
! Number of the cycles for the optimization: |
| 29 |
! XGeomF(NGeomF): Final geometries. |
| 30 |
REAL(KREAL), INTENT(OUT) :: XGeomF(NGeomF) |
| 31 |
|
| 32 |
INTEGER(KINT) :: IdxGeom, I, J, K, Idx, IdxAtom |
| 33 |
REAL(KREAL) :: Rmsd,MRot(3,3), ds, u, v |
| 34 |
REAL(KREAL) :: a_val, d |
| 35 |
|
| 36 |
REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:), DerInt(:) ! (Nat,3) |
| 37 |
REAL(KREAL), ALLOCATABLE :: Xyz_k(:,:) ! (Nat,3) |
| 38 |
REAL(KREAL), ALLOCATABLE :: IntCoord_interpol(:) ! (3*Nat-6) |
| 39 |
REAL(KREAL), ALLOCATABLE :: IntCoord_k(:) ! (3*Nat-6) |
| 40 |
REAL(KREAL), ALLOCATABLE :: XPrimRef(:),XPrim(:) ! NPrim |
| 41 |
|
| 42 |
LOGICAL :: debug, print,printspline |
| 43 |
LOGICAL, EXTERNAL :: valid |
| 44 |
|
| 45 |
INTEGER(KINT) :: NSpline |
| 46 |
CHARACTER(LCHARS) :: FileSpline,TmpChar |
| 47 |
|
| 48 |
|
| 49 |
! We will calculate the length of the path, in MW coordinates... |
| 50 |
! this is done in a stupid way: we interpolate the Baker coordinates values, |
| 51 |
! convert them into cartesian, weight the cartesian |
| 52 |
! and calculate the evolution of the distance ! |
| 53 |
! We have to follow the same procedure for every geometry, |
| 54 |
! so even for the first one, we have to convert it from internal Baker |
| 55 |
! coordinates to cartesian ! |
| 56 |
|
| 57 |
debug=valid("Extrapol_baker")
|
| 58 |
print=valid("printgeom")
|
| 59 |
printspline=(valid("printspline").AND.(dist<=1e30))
|
| 60 |
|
| 61 |
if (debug) WRITE(*,*) "================= Entering Extrapol_baker ====================" |
| 62 |
if (debug) WRITE(*,*) "DBG Extrapol_baker dist=",Dist |
| 63 |
NSpline=int(NMaxPtPath/100) |
| 64 |
!IF (printspline) THEN |
| 65 |
! WRITE(TmpChar,'(I5)') Iopt |
| 66 |
! FileSpline=Trim(adjustL(PathName)) // '_spline.' // AdjustL(TRIM(TmpChar)) |
| 67 |
!OPEN(IOTMP,FILE=FileSpline) |
| 68 |
! END IF |
| 69 |
|
| 70 |
ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),IntCoord_interpol(NCoord),DerInt(NCoord)) |
| 71 |
ALLOCATE(IntCoord_k(NCoord),Xyz_k(Nat,3)) |
| 72 |
ALLOCATE(XPrimRef(NPrim),XPrim(NPrim)) |
| 73 |
|
| 74 |
! XyzGeomI(:,:,:) ! (NGeomI,3,Nat) |
| 75 |
! IntCoordI(:,:) ! (NGeomI,3*Nat-6) |
| 76 |
|
| 77 |
!XyzGeomF(1,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/)) |
| 78 |
XyzGeomF(1,:,:)=XyzGeomI(1,:,:) ! 1st index is geometry-index. |
| 79 |
IntCoordF(1,:)=IntCoordI(1,:) |
| 80 |
|
| 81 |
! We calculate the first derivatives |
| 82 |
u=0.d0 |
| 83 |
DO I=1,NCoord |
| 84 |
! Given the arrays xgeom(1:NGeomI) and IntCoordI(1:NGeomI,Idx) of length |
| 85 |
! NGeomI, which tabulate a function |
| 86 |
! (with the xgeom's in order), and given the array Coef(1:NGeomI,Idx), |
| 87 |
! which is the output from spline, and given a value of u, |
| 88 |
! this routine returns a cubic-spline interpolated value v. |
| 89 |
! and the derivative DerInt(Idx). |
| 90 |
call splintder(u,v,DerInt(I),NGeomI,xgeom(1),IntCoordI(1,I),Coef(1,I)) |
| 91 |
END DO |
| 92 |
IntTangent(1,:)=DerInt |
| 93 |
|
| 94 |
IF (print.AND.(Dist.LE.1e20)) THEN |
| 95 |
WRITE(IOOUT,'(1X,I5)') Nat |
| 96 |
WRITE(IOOUT,*) "# Cartesian Coordinates for geom",1 |
| 97 |
DO I=1,Nat |
| 98 |
If (Renum) THEN |
| 99 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & |
| 100 |
(XyzTmp2(Order(I),J),J=1,3) |
| 101 |
ELSE |
| 102 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), & |
| 103 |
(XyzTmp2(I,J),J=1,3) |
| 104 |
END IF |
| 105 |
END DO |
| 106 |
END IF ! matches IF (print.AND.(Dist.LE.1e20)) THEN |
| 107 |
|
| 108 |
XyzTmp(:,1) = XyzGeomI(1,1,:) ! 1st index is geometry-index. |
| 109 |
XyzTmp(:,2) = XyzGeomI(1,2,:) |
| 110 |
XyzTmp(:,3) = XyzGeomI(1,3,:) |
| 111 |
|
| 112 |
s=0.d0 |
| 113 |
IntCoord_k=IntCoordF(1,:) |
| 114 |
Xyz_k(:,1) = XyzGeomI(1,1,:) ! 1st index is geometry-index. |
| 115 |
Xyz_k(:,2) = XyzGeomI(1,2,:) |
| 116 |
Xyz_k(:,3) = XyzGeomI(1,3,:) |
| 117 |
IdxGeom=1 |
| 118 |
XPrimRef=XPrimitive(1,:) |
| 119 |
XPrimitiveF(1,:)=XPrimitive(1,:) |
| 120 |
DO K=1,NMaxPtPath |
| 121 |
u=real(K)/NMaxPtPath*(NGeomI-1.) |
| 122 |
|
| 123 |
! We generate the interpolated internal coordinates in v. |
| 124 |
! Given the arrays Xgeom(1:NGeomI) (Xgeom(NGeomI): abscissa of all initial geometries) |
| 125 |
! and IntCoordI(1:NGeomI,I) of length NGeomI, which tabulate a function (with the |
| 126 |
! Xgeom's in order), and given the array Coef(1:NGeomI,Idx), which is the output from |
| 127 |
! spline, and given a value of u, this routine returns a cubic-spline interpolated |
| 128 |
! value v and the derivative DerInt(I). |
| 129 |
|
| 130 |
! this loop is to be confirmed: |
| 131 |
! IntCoordI(NGeomI,3*Nat-6) |
| 132 |
DO I=1,NCoord |
| 133 |
call splintder(u,v,DerInt(I),NGeomI,Xgeom(1),IntCoordI(1,I),Coef(1,I)) |
| 134 |
IntCoord_interpol(I)=v |
| 135 |
END DO |
| 136 |
IF(.NOT.FirstTimePathCreate) Then |
| 137 |
WRITE(*,*) "DBG Extrapol_baker Umat_local=UMatF" |
| 138 |
DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart() |
| 139 |
BTransInv_local(I,:) = BTransInvF(IdxGeom,I,:) |
| 140 |
UMat_local(:,I) = UMatF(IdxGeom,:,I) |
| 141 |
END DO |
| 142 |
END IF |
| 143 |
! We convert it into Cartesian coordinates: |
| 144 |
if (debug) WRITE(*,*) "DBG Extrapol_baker, call ConvertBakerInt_car for k=",k |
| 145 |
Call ConvertBakerInternal_cart(IntCoord_k,IntCoord_interpol,Xyz_k(1,1), & |
| 146 |
Xyz_k(1,2),Xyz_k(1,3),XyzTMP2(1,1),XyzTMP2(1,2),XyzTMP2(1,3),XPrim,XPrimRef) |
| 147 |
XPrimRef=Xprim |
| 148 |
IF(.NOT.FirstTimePathCreate) Then |
| 149 |
DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart() |
| 150 |
BTransInvF(IdxGeom,I,:) = BTransInv_local(I,:) |
| 151 |
END DO |
| 152 |
END IF |
| 153 |
|
| 154 |
if (debug) THEN |
| 155 |
WRITE(*,*) "DBG Extrapol_baker, XyzTmp2 before RMSD" |
| 156 |
DO I=1,Nat |
| 157 |
IF (Renum) THEN |
| 158 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & |
| 159 |
(XyzTmp2(Order(I),J),J=1,3) |
| 160 |
ELSE |
| 161 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), & |
| 162 |
(XyzTmp2(I,J),J=1,3) |
| 163 |
END IF |
| 164 |
END DO |
| 165 |
END IF |
| 166 |
|
| 167 |
|
| 168 |
call CalcRmsd(Nat,XyzTmp(1:Nat,1),XyzTmp(1:Nat,2),XyzTmp(1:Nat,3), & |
| 169 |
XyzTmp2(1:Nat,1),XyzTmp2(1:Nat,2),XyzTmp2(1:Nat,3), & |
| 170 |
MRot,rmsd,.TRUE.,.TRUE.) |
| 171 |
|
| 172 |
|
| 173 |
if (debug) THEN |
| 174 |
WRITE(*,*) "DBG Extrapol_baker, XyzTmp2 after RMSD" |
| 175 |
DO I=1,Nat |
| 176 |
IF (Renum) THEN |
| 177 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & |
| 178 |
(XyzTmp2(Order(I),J),J=1,3) |
| 179 |
ELSE |
| 180 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), & |
| 181 |
(XyzTmp2(I,J),J=1,3) |
| 182 |
END IF |
| 183 |
END DO |
| 184 |
END IF |
| 185 |
|
| 186 |
|
| 187 |
IntCoord_k=IntCoord_interpol |
| 188 |
Xyz_k(:,1)=XyzTMP2(:,1) |
| 189 |
Xyz_k(:,2)=XyzTMP2(:,2) |
| 190 |
Xyz_k(:,3)=XyzTMP2(:,3) |
| 191 |
|
| 192 |
ds=0. |
| 193 |
DO I=1,Nat |
| 194 |
DO J=1,3 |
| 195 |
ds=ds+MassAt(I)*(XYZTMp2(I,J)-XYZTmp(I,J))**2 |
| 196 |
XYZTmp(I,J)=XyzTMP2(I,J) |
| 197 |
END DO |
| 198 |
END DO |
| 199 |
|
| 200 |
s=s+sqrt(ds) |
| 201 |
|
| 202 |
IF (s>=dist) THEN |
| 203 |
if (debug) THEN |
| 204 |
WRITE(*,*) "DBG Extrapol_baker s,IdxGeom,dist",s,IdxGeom,dist |
| 205 |
WRITE(*,'(50(1X,F12.8))') IntCoord_interpol |
| 206 |
WRITE(*,*) "DBG Extrapol_baker Angles in deg ?" |
| 207 |
WRITE(*,'(50(1X,F12.8))') IntCoord_interpol*180./pi |
| 208 |
END IF |
| 209 |
s=s-dist |
| 210 |
IdxGeom=IdxGeom+1 |
| 211 |
XprimitiveF(IdxGeom,:)=Xprim(:) |
| 212 |
UMatF(IdxGeom,:,:)=UMat_local(:,:) |
| 213 |
SGeom(IdxGeom)=s+IdxGeom*dist !SGeom(NGeomF) |
| 214 |
XgeomF(IdxGeom)=u |
| 215 |
XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/)) |
| 216 |
|
| 217 |
! IntCoordF(NGeomF,NCoord): Final Internal coordinates for number of final |
| 218 |
! geometries. NCoord is the number of coordinates (NCoord) of each geometry. |
| 219 |
IntCoordF(IdxGeom,:)=IntCoord_interpol(:) |
| 220 |
IntTangent(IdxGeom,:)=DerInt |
| 221 |
|
| 222 |
IF (print) THEN |
| 223 |
WRITE(IOOUT,'(1X,I5)') Nat |
| 224 |
WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K |
| 225 |
! PFL 17/July/2006: only if we have more than 4 atoms. |
| 226 |
IF (Nat.GE.4) THEN |
| 227 |
Call CalcRmsd(Nat,x0,y0,z0, & |
| 228 |
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & |
| 229 |
MRot,rmsd,.TRUE.,.TRUE.) |
| 230 |
END IF |
| 231 |
|
| 232 |
DO I=1,Nat |
| 233 |
IF (Renum) THEN |
| 234 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & |
| 235 |
(XyzTmp2(Order(I),J),J=1,3) |
| 236 |
ELSE |
| 237 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), & |
| 238 |
(XyzTmp2(I,J),J=1,3) |
| 239 |
END IF |
| 240 |
END DO |
| 241 |
END IF !matches IF (print) THEN |
| 242 |
END IF ! matches IF (s>=dist) THEN |
| 243 |
END DO ! matches DO K=1,NMaxPtPath |
| 244 |
|
| 245 |
|
| 246 |
if (s>=0.9*dist) THEN |
| 247 |
s=s-dist |
| 248 |
IdxGeom=IdxGeom+1 |
| 249 |
SGeom(IdxGeom)=s+IdxGeom*dist |
| 250 |
XgeomF(IdxGeom)=min(u,NGeomI-1.d0) |
| 251 |
XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/)) |
| 252 |
! XyzGeomF(IdxGeom,:,:)=XyzTmp2(:,:) |
| 253 |
|
| 254 |
IntCoordF(IdxGeom,:)=IntCoord_interpol(:) |
| 255 |
XprimitiveF(IdxGeom,:)=Xprim(:) |
| 256 |
UMatF(IdxGeom,:,:)=UMat_local(:,:) |
| 257 |
IntTangent(IdxGeom,:)=DerInt |
| 258 |
|
| 259 |
if (print) THEN |
| 260 |
WRITE(IOOUT,'(1X,I5)') Nat |
| 261 |
WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K |
| 262 |
! PFL 17/July/2006: only if we have more than 4 atoms. |
| 263 |
IF (Nat.GE.4) THEN |
| 264 |
Call CalcRmsd(Nat,x0,y0,z0, & |
| 265 |
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & |
| 266 |
MRot,rmsd,.TRUE.,.TRUE.) |
| 267 |
END IF |
| 268 |
|
| 269 |
DO I=1,Nat |
| 270 |
IF (Renum) THEN |
| 271 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & |
| 272 |
(XyzTmp2(Order(I),J),J=1,3) |
| 273 |
ELSE |
| 274 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), & |
| 275 |
(XyzTmp2(I,J),J=1,3) |
| 276 |
END IF |
| 277 |
END DO |
| 278 |
END IF ! matches if (print) THEN |
| 279 |
END IF ! matches if (s>=0.9*dist) THEN |
| 280 |
|
| 281 |
if (debug) WRITE(*,*) 's final =',s |
| 282 |
if (debug) THEN |
| 283 |
WRITE(*,*) "XPrimitiveF" |
| 284 |
DO I=1,NGeomF |
| 285 |
WRITE(*,'(1X,I5," : ",50(1X,F10.6))') I,XPrimitiveF(I,:) |
| 286 |
END DO |
| 287 |
END IF |
| 288 |
|
| 289 |
DEALLOCATE(XyzTmp,XyzTmp2,IntCoord_interpol,IntCoord_k,Xyz_k) |
| 290 |
|
| 291 |
if (printspline) CLOSE(IOTMP) |
| 292 |
if (debug) WRITE(*,*) "================= Extrapol_baker Over =====================" |
| 293 |
|
| 294 |
END SUBROUTINE EXTRAPOL_BAKER |