Révision 8
| src/Egrad_baker.f90 (revision 8) | ||
|---|---|---|
| 1 |
subroutine Egrad_baker(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom,GeomOld,GeomCart_old) |
|
| 2 |
|
|
| 3 |
! This routines calculates the energy E and the gradient Grad of |
|
| 4 |
! a molecule with Geometry Geom (may be in internal coordinates), |
|
| 5 |
! using for now, either Gaussian or Ext, more general later. |
|
| 6 |
|
|
| 7 |
use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, & |
|
| 8 |
prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, & |
|
| 9 |
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & |
|
| 10 |
, BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & |
|
| 11 |
, Order,OrderInv, XPrimitiveF |
|
| 12 |
|
|
| 13 |
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord) |
|
| 14 |
! allocated in Path.f90 |
|
| 15 |
|
|
| 16 |
use Io_module |
|
| 17 |
IMPLICIT NONE |
|
| 18 |
|
|
| 19 |
! Energy (calculated if F300K=.F., else estimated) |
|
| 20 |
REAL(KREAL), INTENT (OUT) :: E |
|
| 21 |
! NCoord: Number of the degrees of freedom |
|
| 22 |
! IGeom: index of the geometry. |
|
| 23 |
INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt |
|
| 24 |
! Geometry at which gradient is calculated (cf Factual also): |
|
| 25 |
REAL(KREAL), INTENT (INOUT) :: Geom(NCoord) |
|
| 26 |
REAL(KREAL), INTENT (INOUT) :: GeomOld(NCoord) |
|
| 27 |
! Gradient calculated at Geom geometry: |
|
| 28 |
REAL(KREAL), INTENT (OUT) :: Grad(NCoord) |
|
| 29 |
! Cartesian geometry corresponding to (Internal Geometry) Geom: |
|
| 30 |
REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3) |
|
| 31 |
REAL(KREAL), INTENT (IN) :: GeomCart_old(Nat,3) |
|
| 32 |
LOGICAL,INTENT (IN) :: Flag_Opt_Geom |
|
| 33 |
|
|
| 34 |
! ====================================================================== |
|
| 35 |
|
|
| 36 |
CHARACTER(LCHARS) :: LINE |
|
| 37 |
LOGICAL :: debug |
|
| 38 |
LOGICAL, SAVE :: first=.true. |
|
| 39 |
LOGICAL, ALLOCATABLE :: Done(:) |
|
| 40 |
REAL(KREAL), SAVE :: Eold=1.e6 |
|
| 41 |
REAL(KREAL), ALLOCATABLE :: GradTmp(:) |
|
| 42 |
REAL(KREAL), ALLOCATABLE :: GradCart(:) |
|
| 43 |
REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) |
|
| 44 |
REAL(KREAL) :: d, a_val, Pi |
|
| 45 |
|
|
| 46 |
INTEGER(KINT) :: iat, jat, kat, i, j, n3at,IBeg,Idx |
|
| 47 |
INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11 |
|
| 48 |
|
|
| 49 |
CHARACTER(LCHARS) :: ListName, CH32SVAR1 |
|
| 50 |
CHARACTER(VLCHARS), SAVE :: RstrtCopy, RunCommand |
|
| 51 |
LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False. |
|
| 52 |
LOGICAL :: FTmp |
|
| 53 |
INTEGER(kint) :: NbLists,LastIt,Nat4,NStep, Firstit |
|
| 54 |
INTEGER(KINT), PARAMETER :: NbExtName=4 |
|
| 55 |
|
|
| 56 |
REAL(KREAL), ALLOCATABLE :: XPrimRef(:), XPrim(:) ! NPrim |
|
| 57 |
|
|
| 58 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 59 |
! |
|
| 60 |
!!!!!!!! PFL Debug |
|
| 61 |
! |
|
| 62 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 63 |
REAL(KREAL), ALLOCATABLE :: GeomTmp(:,:) !3,Nat |
|
| 64 |
REAL(KREAL), ALLOCATABLE :: V(:,:) ! (NCoord,NCoord) |
|
| 65 |
REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:) ! xyzatm(3,nat) |
|
| 66 |
REAL(KREAL) :: tol |
|
| 67 |
LOGICAL :: FAIL , FALOBBT,FALOBPrimt,FAloBBTInv |
|
| 68 |
LOGICAL :: DebugPFL |
|
| 69 |
|
|
| 70 |
! ====================================================================== |
|
| 71 |
LOGICAL, EXTERNAL :: valid |
|
| 72 |
! ====================================================================== |
|
| 73 |
|
|
| 74 |
INTEGER(KINT) :: OptGeom |
|
| 75 |
Namelist/path/OptGeom |
|
| 76 |
|
|
| 77 |
Pi=dacos(-1.0d0) |
|
| 78 |
n3at=3*nat |
|
| 79 |
|
|
| 80 |
debug=valid('EGRAD_baker')
|
|
| 81 |
debugPFL=valid('bakerPFL')
|
|
| 82 |
if (debug) WRITE(*,*) '================ Entering Egrad_baker =================' |
|
| 83 |
if (debug) Print *, 'EGrad.f90, L73, IOpt=', IOpt, ' Flag_Opt_Geom=',Flag_Opt_Geom |
|
| 84 |
|
|
| 85 |
Grad=0. |
|
| 86 |
E=0. |
|
| 87 |
|
|
| 88 |
ALLOCATE(GradCart(3*Nat)) |
|
| 89 |
ALLOCATE(x(Nat),y(Nat),z(Nat)) |
|
| 90 |
ALLOCATE(XPrimRef(NPrim),XPrim(NPrim)) |
|
| 91 |
|
|
| 92 |
SELECT CASE (COORD) |
|
| 93 |
CASE ('ZMAT')
|
|
| 94 |
WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad.f90, instead of Egrad_baker.f90. Stop" |
|
| 95 |
CASE ('BAKER')
|
|
| 96 |
XPrimRef=XPrimitiveF(IGeom,:) |
|
| 97 |
IF (Flag_Opt_Geom) THEN |
|
| 98 |
IF (IOpt .EQ. 1) THEN |
|
| 99 |
GeomCart(:,1) = XyzGeomI(IGeom,1,:) |
|
| 100 |
GeomCart(:,2) = XyzGeomI(IGeom,2,:) |
|
| 101 |
GeomCart(:,3) = XyzGeomI(IGeom,3,:) |
|
| 102 |
ELSE |
|
| 103 |
if (debug) WRITE(*,*) 'Before Call ConvertBakerInternal_cart; L82, Egrad.f90' |
|
| 104 |
WRITE(*,*) 'GeomOld=' |
|
| 105 |
WRITE(*,'(12(1X,F6.3))') GeomOld |
|
| 106 |
WRITE(*,*) 'Geom=' |
|
| 107 |
WRITE(*,'(12(1X,F6.3))') Geom |
|
| 108 |
WRITE(*,*) 'GeomCart_old=' |
|
| 109 |
WRITE(*,'(20(1X,F6.3))') GeomCart_old |
|
| 110 |
Call ConvertBakerInternal_cart(GeomOld,Geom,GeomCart_old(:,1), & |
|
| 111 |
GeomCart_old(:,2),GeomCart_old(:,3),x,y,z,XPrim,XPrimRef) |
|
| 112 |
GeomCart(:,1) = x(:) |
|
| 113 |
GeomCart(:,2) = y(:) |
|
| 114 |
GeomCart(:,3) = z(:) |
|
| 115 |
END IF |
|
| 116 |
ELSE |
|
| 117 |
IF (IOpt .EQ. 0) THEN |
|
| 118 |
GeomCart(:,1) = XyzGeomF(IGeom,1,:) |
|
| 119 |
GeomCart(:,2) = XyzGeomF(IGeom,2,:) |
|
| 120 |
GeomCart(:,3) = XyzGeomF(IGeom,3,:) |
|
| 121 |
ELSE |
|
| 122 |
! IntCoordF(NGeomF,NCoord), XyzGeomF(NGeomF,3,Nat) |
|
| 123 |
! Geom has to be converted into x,y,z |
|
| 124 |
! GeomOld_all and XyzGeomF are old geometry in internal and cartesian coords resp. |
|
| 125 |
DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart() |
|
| 126 |
BTransInv_local(I,:) = BTransInvF(IGeom,I,:) |
|
| 127 |
UMat_local(:,I) = UMatF(IGeom,:,I) |
|
| 128 |
END DO |
|
| 129 |
if (debug) Print *, 'I am before Call ConvertBakerInternal_cart; L108, Egrad.f90' |
|
| 130 |
WRITE(*,*) 'GeomOld=' |
|
| 131 |
WRITE(*,'(12(1X,F6.3))') GeomOld |
|
| 132 |
WRITE(*,*) 'Geom=' |
|
| 133 |
WRITE(*,'(12(1X,F6.3))') Geom |
|
| 134 |
WRITE(*,*) 'DBG Egrad_baker L129 GeomCart_old=' |
|
| 135 |
DO I=1,Nat |
|
| 136 |
WRITE(*,'(3(1X,F12.8))') XyzGeomF(IGeom,:,I) |
|
| 137 |
END DO |
|
| 138 |
|
|
| 139 |
Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),Geom,XyzGeomF(IGeom,1,:), & |
|
| 140 |
XyzGeomF(IGeom,2,:),XyzGeomF(IGeom,3,:),x,y,z,XPrim,XPrimRef) |
|
| 141 |
if (debug) Print *, 'I am after Call ConvertBakerInternal_cart; L111, Egrad.f90' |
|
| 142 |
DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart() |
|
| 143 |
BTransInvF(IGeom,I,:) = BTransInv_local(I,:) |
|
| 144 |
END DO |
|
| 145 |
! GeomCart is actually XyzCoord and has INTENT(OUT) in this subroutine. |
|
| 146 |
GeomCart(:,1) = x(:) |
|
| 147 |
GeomCart(:,2) = y(:) |
|
| 148 |
GeomCart(:,3) = z(:) |
|
| 149 |
END IF ! matches IF (IOpt .EQ. 0) THEN |
|
| 150 |
END IF |
|
| 151 |
CASE ('MIXED')
|
|
| 152 |
WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad.f90, instead of Egrad_baker.f90. Stop" |
|
| 153 |
CASE ('CART')
|
|
| 154 |
WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad.f90, instead of Egrad_baker.f90. Stop" |
|
| 155 |
CASE ('HYBRID')
|
|
| 156 |
WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad.f90, instead of Egrad_baker.f90. Stop" |
|
| 157 |
CASE DEFAULT |
|
| 158 |
WRITE (*,*) "Coord=",TRIM(COORD)," not recognized. Stop" |
|
| 159 |
STOP |
|
| 160 |
END SELECT |
|
| 161 |
!WRITE(*,*) Nat |
|
| 162 |
!WRITE(*,*) 'GeomCart in Egrad_baker' |
|
| 163 |
!DO I=1,Nat |
|
| 164 |
! WRITE(*,'(1X,A,3(1X,F12.8))') AtName(I),GeomCart(I,1:3) |
|
| 165 |
!END DO |
|
| 166 |
|
|
| 167 |
SELECT CASE (Prog) |
|
| 168 |
CASE ('GAUSSIAN')
|
|
| 169 |
! we call the Gaussian routine. |
|
| 170 |
! it will return the energy and the gradient in cartesian coordinates. |
|
| 171 |
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT) |
|
| 172 |
Call egrad_gaussian(E,GeomCart,GradCart) |
|
| 173 |
CASE ('MOPAC')
|
|
| 174 |
! we call the Mopac routine. |
|
| 175 |
! it will return the energy and the gradient in cartesian coordinates. |
|
| 176 |
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT) |
|
| 177 |
Call egrad_mopac(E,GeomCart,GradCart) |
|
| 178 |
CASE ('VASP')
|
|
| 179 |
Call egrad_vasp(E,Geomcart,GradCart) |
|
| 180 |
CASE ('TURBOMOLE')
|
|
| 181 |
Call egrad_turbomole(E,Geomcart,GradCart) |
|
| 182 |
CASE ('EXT')
|
|
| 183 |
Call egrad_ext(E,Geomcart,GradCart) |
|
| 184 |
CASE ('TEST')
|
|
| 185 |
Call egrad_test(Nat,E,Geomcart,GradCart) |
|
| 186 |
CASE ('CHAMFRE')
|
|
| 187 |
Call egrad_chamfre(Nat,E,Geomcart,GradCart) |
|
| 188 |
END SELECT |
|
| 189 |
if (debug) THEN |
|
| 190 |
WRITE(*,*) 'DBG EGRAD, GradCart read' |
|
| 191 |
DO Iat=1,Nat |
|
| 192 |
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GradCart(3*Iat-2:3*Iat) |
|
| 193 |
END DO |
|
| 194 |
END IF |
|
| 195 |
|
|
| 196 |
! If (PROG=="VASP") STOP |
|
| 197 |
|
|
| 198 |
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid |
|
| 199 |
! but we have to convert it into Zmat if COORD=Zmat |
|
| 200 |
SELECT CASE (COORD) |
|
| 201 |
CASE ("ZMAT")
|
|
| 202 |
WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad.f90, instead of Egrad_baker.f90. Stop" |
|
| 203 |
CASE ('BAKER')
|
|
| 204 |
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid |
|
| 205 |
! but we have to convert it into internal coordinates if COORD=Baker |
|
| 206 |
!!!!!!!!!!!!!!!!!!!! |
|
| 207 |
! |
|
| 208 |
! PFL Debug |
|
| 209 |
! |
|
| 210 |
!!!!!!!!!!!!!!!!!!!! |
|
| 211 |
if (debugPFL) THEN |
|
| 212 |
WRITE(*,*) "DBG PFL Egrad_Baker... computing dzdc" |
|
| 213 |
if (.not.ALLOCATED(indzmat)) THEN |
|
| 214 |
ALLOCATE(indzmat(Nat,5)) |
|
| 215 |
IndZmat=0 |
|
| 216 |
Indzmat(1,1)=1 |
|
| 217 |
IndZmat(2,1)=2 |
|
| 218 |
IndZmat(2,2)=1 |
|
| 219 |
IndZmat(3,1)=3 |
|
| 220 |
IndZmat(3,2)=2 |
|
| 221 |
IndZmat(3,3)=1 |
|
| 222 |
IndZmat(4,1)=4 |
|
| 223 |
IndZmat(4,2)=3 |
|
| 224 |
IndZmat(4,3)=2 |
|
| 225 |
IndZmat(4,4)=1 |
|
| 226 |
END IF |
|
| 227 |
IF (.NOT.ALLOCATED(DzDc)) THEN |
|
| 228 |
ALLOCATE(DzDc(3,nat,3,Nat)) |
|
| 229 |
END IF |
|
| 230 |
if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat)) |
|
| 231 |
DO I=1,Nat |
|
| 232 |
GeomTmp(1,I)=GeomCart(OrderInv(I),1) |
|
| 233 |
GeomTmp(2,I)=GeomCart(OrderInv(i),2) |
|
| 234 |
GeomTmp(3,I)=GeomCart(OrderInv(i),3) |
|
| 235 |
END DO |
|
| 236 |
|
|
| 237 |
DzDc=0.d0 |
|
| 238 |
CALL CalcBMat_int (nat, GeomTmp, indzmat, dzdc, massat,atome) |
|
| 239 |
|
|
| 240 |
END IF ! if (debugPFL) |
|
| 241 |
!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 242 |
! Debugging PFL |
|
| 243 |
!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 244 |
|
|
| 245 |
WRITE(*,*) "BTransInv_local Trans that is used originally" |
|
| 246 |
DO J=1,3*Nat |
|
| 247 |
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) |
|
| 248 |
END DO |
|
| 249 |
|
|
| 250 |
WRITE(*,*) "Calculating actual values using GeomCart" |
|
| 251 |
if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat)) |
|
| 252 |
GeomTmp(1,:)=GeomCart(1:Nat,1) |
|
| 253 |
GeomTmp(2,:)=GeomCart(1:Nat,2) |
|
| 254 |
GeomTmp(3,:)=GeomCart(1:Nat,3) |
|
| 255 |
|
|
| 256 |
! Computing B^prim: |
|
| 257 |
FAloBBT=.NOT.ALLOCATED(BBT) |
|
| 258 |
FAloBBTInv=.NOT.ALLOCATED(BBT_inv) |
|
| 259 |
FAloBPrimT=.NOT.ALLOCATED(BprimT) |
|
| 260 |
if (FAloBPrimT) ALLOCATE(Bprimt(3*Nat,NPrim)) |
|
| 261 |
if (FAloBBT) ALLOCATE(BBt(NCoord,NCoord)) |
|
| 262 |
if (FAloBBTInv) ALLOCATE(BBt_inv(NCoord,NCoord)) |
|
| 263 |
BprimT=0.d0 |
|
| 264 |
ScanCoord=>Coordinate |
|
| 265 |
I=0 |
|
| 266 |
DO WHILE (Associated(ScanCoord%next)) |
|
| 267 |
I=I+1 |
|
| 268 |
SELECT CASE (ScanCoord%Type) |
|
| 269 |
CASE ('BOND')
|
|
| 270 |
CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,GeomTmp,BprimT(1,I)) |
|
| 271 |
CASE ('ANGLE')
|
|
| 272 |
CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,GeomTmp,BprimT(1,I)) |
|
| 273 |
CASE ('DIHEDRAL')
|
|
| 274 |
CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,GeomTmp,BprimT(1,I)) |
|
| 275 |
END SELECT |
|
| 276 |
ScanCoord => ScanCoord%next |
|
| 277 |
END DO |
|
| 278 |
|
|
| 279 |
if (debug) THEN |
|
| 280 |
WRITE(*,*) "Bprim " |
|
| 281 |
DO J=1,3*Nat |
|
| 282 |
WRITE(*,'(50(1X,F12.6))') BprimT(J,:) |
|
| 283 |
END DO |
|
| 284 |
END IF |
|
| 285 |
|
|
| 286 |
BMat_BakerT = 0.d0 |
|
| 287 |
DO I=1,NCoord |
|
| 288 |
DO J=1,NPrim |
|
| 289 |
!BprimT is transpose of B^prim. |
|
| 290 |
!B = UMat^T B^prim, B^T = (B^prim)^T UMat |
|
| 291 |
BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I) |
|
| 292 |
END DO |
|
| 293 |
END DO |
|
| 294 |
|
|
| 295 |
BBT = 0.d0 |
|
| 296 |
DO I=1,NCoord |
|
| 297 |
DO J=1,3*Nat ! BBT(:,I) forms BB^T |
|
| 298 |
BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I) |
|
| 299 |
END DO |
|
| 300 |
END DO |
|
| 301 |
|
|
| 302 |
BBT_inv = 0.d0 |
|
| 303 |
|
|
| 304 |
Call GenInv(NCoord,BBT,BBT_inv,NCoord) |
|
| 305 |
|
|
| 306 |
! ALLOCATE(V(NCoord,NCoord)) |
|
| 307 |
! tol = 1e-10 |
|
| 308 |
! ! BBT is destroyed by GINVSE |
|
| 309 |
! Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT) |
|
| 310 |
! DEALLOCATE(V) |
|
| 311 |
! IF(Fail) Then |
|
| 312 |
! WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90" |
|
| 313 |
! STOP |
|
| 314 |
! END IF |
|
| 315 |
|
|
| 316 |
!Print *, 'BBT_inv=' |
|
| 317 |
DO J=1,NCoord |
|
| 318 |
! WRITE(*,'(50(1X,F10.2))') BBT_inv(:,J) |
|
| 319 |
! Print *, BBT_inv(:,J) |
|
| 320 |
END DO |
|
| 321 |
!Print *, 'End of BBT_inv' |
|
| 322 |
|
|
| 323 |
! Calculation of (B^T)^-1 = (BB^T)^-1B: |
|
| 324 |
BTransInv_local = 0.d0 |
|
| 325 |
DO I=1,3*Nat |
|
| 326 |
DO J=1,NCoord |
|
| 327 |
BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J) |
|
| 328 |
END DO |
|
| 329 |
END DO |
|
| 330 |
|
|
| 331 |
Print *, 'BMat_BakerT=' |
|
| 332 |
DO J=1,3*Nat |
|
| 333 |
WRITE(*,'(50(1X,F12.8))') BMat_BakerT(:,J) |
|
| 334 |
END DO |
|
| 335 |
|
|
| 336 |
WRITE(*,*) "BTransInv_local Trans corresponding to GeomCart" |
|
| 337 |
DO J=1,3*Nat |
|
| 338 |
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) |
|
| 339 |
END DO |
|
| 340 |
|
|
| 341 |
if (FAloBPrimT) DEALLOCATE(Bprimt) |
|
| 342 |
if (FAloBBT) DEALLOCATE(BBt) |
|
| 343 |
if (FAloBBTInv) DEALLOCATE(BBt_inv) |
|
| 344 |
|
|
| 345 |
Grad=0.d0 |
|
| 346 |
DO I=1, 3*Nat |
|
| 347 |
! Here Grad is gradient in Baker coordinates. GradCart(3*Nat) |
|
| 348 |
Grad(:) = Grad(:) + BTransInv_local(:,I)*GradCart(I) !BTransInv_local = BTransInv in Opt_Geom.f90 |
|
| 349 |
END DO |
|
| 350 |
!WRITE(*,*) 'In Egrad.f90, L263, Gradient:' |
|
| 351 |
!WRITE(*,'(12(1X,F8.4))') Grad(1:NCoord) |
|
| 352 |
CASE ('MIXED')
|
|
| 353 |
WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad.f90, instead of Egrad_baker.f90. Stop" |
|
| 354 |
CASE ("CART","HYBRID")
|
|
| 355 |
WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad.f90, instead of Egrad_baker.f90. Stop" |
|
| 356 |
CASE DEFAULT |
|
| 357 |
WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in Egradients. STOP" |
|
| 358 |
STOP |
|
| 359 |
END SELECT |
|
| 360 |
|
|
| 361 |
DEALLOCATE(GradCart) |
|
| 362 |
DEALLOCATE(x,y,z) |
|
| 363 |
IF (Debug) WRITE(*,*) 'DBG EGRAD_baker grad',grad(1:NCoord) |
|
| 364 |
|
|
| 365 |
if (debug) WRITE(*,*) '================ Egrad_baker Over ====================' |
|
| 366 |
|
|
| 367 |
!999 CONTINUE |
|
| 368 |
!if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!' |
|
| 369 |
!STOP |
|
| 370 |
! ====================================================================== |
|
| 371 |
end subroutine Egrad_baker |
|
| 372 |
|
|
| src/Egrad.f90 (revision 8) | ||
|---|---|---|
| 1 |
subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart) |
|
| 1 |
subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old)
|
|
| 2 | 2 |
|
| 3 | 3 |
! This routines calculates the energy E and the gradient Grad of |
| 4 | 4 |
! a molecule with Geometry Geom (may be in internal coordinates), |
| 5 |
! using for now, either gaussian or Ext, more general later.
|
|
| 5 |
! using for now, either Gaussian or Ext, more general later.
|
|
| 6 | 6 |
|
| 7 | 7 |
use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, & |
| 8 |
prog,NCart,XyzGeomF,IntCoordF,XyzGeomI, renum, OrderInv |
|
| 9 |
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord) |
|
| 10 |
! allocated in Path.f90 |
|
| 8 |
prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, renum, & |
|
| 9 |
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & |
|
| 10 |
, BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & |
|
| 11 |
, Order,OrderInv, XPrimitiveF |
|
| 12 |
|
|
| 13 |
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord) |
|
| 14 |
! allocated in Path.f90 |
|
| 15 |
|
|
| 11 | 16 |
use Io_module |
| 12 | 17 |
IMPLICIT NONE |
| 13 | 18 |
|
| ... | ... | |
| 22 | 27 |
REAL(KREAL), INTENT (OUT) :: Grad(NCoord) |
| 23 | 28 |
! Cartesian geometry corresponding to (Internal Geometry) Geom: |
| 24 | 29 |
REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3) |
| 30 |
!!! Optional, just for geometry optimization with Baker coordinates |
|
| 31 |
REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3) |
|
| 32 |
REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord) |
|
| 33 |
! FOptGeom is a flag indicating if we are doing a geom optimization |
|
| 34 |
! it can be omitted so that we use a local flag: Flag_Opt_Geom |
|
| 35 |
LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom |
|
| 36 |
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F |
|
| 37 |
LOGICAL :: Flag_Opt_Geom |
|
| 25 | 38 |
|
| 26 | 39 |
! ====================================================================== |
| 27 | 40 |
|
| 28 |
CHARACTER(LCHARS) :: LINE |
|
| 29 | 41 |
LOGICAL :: debug |
| 30 |
LOGICAL, SAVE :: first=.true. |
|
| 31 |
REAL(KREAL), SAVE :: Eold=1.e6 |
|
| 32 | 42 |
REAL(KREAL), ALLOCATABLE :: GradTmp(:) |
| 33 | 43 |
REAL(KREAL), ALLOCATABLE :: GradCart(:) |
| 34 | 44 |
REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) |
| 35 |
REAL(KREAL) :: d, a_val, Pi
|
|
| 45 |
REAL(KREAL) :: Pi |
|
| 36 | 46 |
|
| 37 |
INTEGER(KINT) :: iat,jat,kat,i,j,IBeg,Idx
|
|
| 47 |
INTEGER(KINT) :: iat, i, j, IBeg,Idx
|
|
| 38 | 48 |
INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11 |
| 39 | 49 |
|
| 40 |
CHARACTER(LCHARS) :: ListName, CH32SVAR1 |
|
| 41 |
CHARACTER(VLCHARS), SAVE :: RstrtCopy, RunCommand |
|
| 42 |
LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False. |
|
| 43 |
LOGICAL :: FTmp |
|
| 44 |
INTEGER(kint) :: NbLists,NStep |
|
| 45 | 50 |
INTEGER(KINT), PARAMETER :: NbExtName=4 |
| 46 | 51 |
|
| 52 |
REAL(KREAL), ALLOCATABLE :: XPrimRef(:), XPrim(:) ! NPrim |
|
| 47 | 53 |
|
| 54 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 55 |
! |
|
| 56 |
!!!!!!!! PFL Debug |
|
| 57 |
! |
|
| 58 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 59 |
REAL(KREAL), ALLOCATABLE :: GeomTmp(:,:) !3,Nat |
|
| 60 |
LOGICAL :: FALOBBT,FALOBPrimt,FAloBBTInv |
|
| 61 |
LOGICAL :: DebugPFL |
|
| 62 |
|
|
| 48 | 63 |
! ====================================================================== |
| 49 | 64 |
INTERFACE |
| 50 | 65 |
function valid(string) result (isValid) |
| ... | ... | |
| 53 | 68 |
END function VALID |
| 54 | 69 |
END INTERFACE |
| 55 | 70 |
! ====================================================================== |
| 56 |
|
|
| 57 |
INTEGER(KINT) :: OptGeom |
|
| 58 |
Namelist/path/OptGeom |
|
| 59 |
|
|
| 71 |
|
|
| 60 | 72 |
Pi=dacos(-1.0d0) |
| 61 | 73 |
|
| 74 |
if (present(FOptGeom)) THEN |
|
| 75 |
Flag_Opt_Geom=FOptGeom |
|
| 76 |
ELSE |
|
| 77 |
Flag_Opt_Geom=.FALSE. |
|
| 78 |
END IF |
|
| 79 |
|
|
| 62 | 80 |
debug=valid('EGRAD')
|
| 63 |
if (debug) Call header("Entering Egrad")
|
|
| 81 |
debugPFL=valid('bakerPFL')
|
|
| 82 |
if (debug) Call Header("Entering Egrad")
|
|
| 64 | 83 |
|
| 65 |
Print *, 'EGrad.f90, L73, IOpt=', IOpt
|
|
| 84 |
if (debug) Print *, 'EGrad.f90, L73, IOpt=', IOpt, ' Flag_Opt_Geom=',Flag_Opt_Geom
|
|
| 66 | 85 |
|
| 67 | 86 |
Grad=0. |
| 68 | 87 |
E=0. |
| 69 | 88 |
|
| 89 |
|
|
| 70 | 90 |
ALLOCATE(GradCart(3*Nat)) |
| 71 | 91 |
ALLOCATE(x(Nat),y(Nat),z(Nat)) |
| 92 |
ALLOCATE(XPrimRef(NPrim),XPrim(NPrim)) |
|
| 93 |
|
|
| 72 | 94 |
SELECT CASE (COORD) |
| 73 | 95 |
CASE ('ZMAT')
|
| 74 | 96 |
! In order to avoid problem with numbering and linear angles/diehedral, we convert the |
| ... | ... | |
| 76 | 98 |
! A remplacer par Int2Cart :-) |
| 77 | 99 |
Call Int2Cart(nat,indzmat,Geom,GeomCart) |
| 78 | 100 |
CASE ('BAKER')
|
| 79 |
WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad_baker.f90, instead of Egrad.f90. Stop" |
|
| 101 |
XPrimRef=XPrimitiveF(IGeom,:) |
|
| 102 |
IF (Flag_Opt_Geom) THEN |
|
| 103 |
IF (IOpt .EQ. 1) THEN |
|
| 104 |
GeomCart(:,1) = XyzGeomI(IGeom,1,:) |
|
| 105 |
GeomCart(:,2) = XyzGeomI(IGeom,2,:) |
|
| 106 |
GeomCart(:,3) = XyzGeomI(IGeom,3,:) |
|
| 107 |
ELSE |
|
| 108 |
if (debug) WRITE(*,*) 'Before Call ConvertBakerInternal_cart; L125, Egrad.f90' |
|
| 109 |
if (.NOT.present(GeomOld)) THEN |
|
| 110 |
WRITE(*,*) "ERROR: GeomOld not passed as an argument" |
|
| 111 |
STOP |
|
| 112 |
END IF |
|
| 113 |
WRITE(*,*) 'GeomOld=' |
|
| 114 |
WRITE(*,'(12(1X,F6.3))') GeomOld |
|
| 115 |
WRITE(*,*) 'Geom=' |
|
| 116 |
WRITE(*,'(12(1X,F6.3))') Geom |
|
| 117 |
if (.NOT.present(GeomCart_Old)) THEN |
|
| 118 |
WRITE(*,*) "ERROR: GeomCart_Old not passed as an argument" |
|
| 119 |
STOP |
|
| 120 |
END IF |
|
| 121 |
|
|
| 122 |
WRITE(*,*) 'GeomCart_old=' |
|
| 123 |
WRITE(*,'(20(1X,F6.3))') GeomCart_old(1:Nat,1:3) |
|
| 124 |
Call ConvertBakerInternal_cart(GeomOld,Geom,GeomCart_old(:,1), & |
|
| 125 |
GeomCart_old(:,2),GeomCart_old(:,3),x,y,z,XPrim,XPrimRef) |
|
| 126 |
GeomCart(:,1) = x(:) |
|
| 127 |
GeomCart(:,2) = y(:) |
|
| 128 |
GeomCart(:,3) = z(:) |
|
| 129 |
END IF |
|
| 130 |
ELSE |
|
| 131 |
IF (IOpt .EQ. 0) THEN |
|
| 132 |
GeomCart(:,1) = XyzGeomF(IGeom,1,:) |
|
| 133 |
GeomCart(:,2) = XyzGeomF(IGeom,2,:) |
|
| 134 |
GeomCart(:,3) = XyzGeomF(IGeom,3,:) |
|
| 135 |
ELSE |
|
| 136 |
! IntCoordF(NGeomF,NCoord), XyzGeomF(NGeomF,3,Nat) |
|
| 137 |
! Geom has to be converted into x,y,z |
|
| 138 |
! GeomOld_all and XyzGeomF are old geometry in internal and cartesian coords resp. |
|
| 139 |
DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart() |
|
| 140 |
BTransInv_local(I,:) = BTransInvF(IGeom,I,:) |
|
| 141 |
UMat_local(:,I) = UMatF(IGeom,:,I) |
|
| 142 |
END DO |
|
| 143 |
if (debug) Print *, 'I am before Call ConvertBakerInternal_cart; L160, Egrad.f90' |
|
| 144 |
WRITE(*,*) 'GeomOld=' |
|
| 145 |
WRITE(*,'(12(1X,F6.3))') GeomOld |
|
| 146 |
WRITE(*,*) 'Geom=' |
|
| 147 |
WRITE(*,'(12(1X,F6.3))') Geom |
|
| 148 |
WRITE(*,*) 'DBG Egrad L165 GeomCart_old=' |
|
| 149 |
DO I=1,Nat |
|
| 150 |
WRITE(*,'(3(1X,F12.8))') XyzGeomF(IGeom,:,I) |
|
| 151 |
END DO |
|
| 152 |
|
|
| 153 |
Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),Geom,XyzGeomF(IGeom,1,:), & |
|
| 154 |
XyzGeomF(IGeom,2,:),XyzGeomF(IGeom,3,:),x,y,z,XPrim,XPrimRef) |
|
| 155 |
if (debug) Print *, 'I am after Call ConvertBakerInternal_cart; L111, Egrad.f90' |
|
| 156 |
DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart() |
|
| 157 |
BTransInvF(IGeom,I,:) = BTransInv_local(I,:) |
|
| 158 |
END DO |
|
| 159 |
! GeomCart is actually XyzCoord and has INTENT(OUT) in this subroutine. |
|
| 160 |
GeomCart(:,1) = x(:) |
|
| 161 |
GeomCart(:,2) = y(:) |
|
| 162 |
GeomCart(:,3) = z(:) |
|
| 163 |
END IF ! matches IF (IOpt .EQ. 0) THEN |
|
| 164 |
END IF |
|
| 80 | 165 |
CASE ('MIXED')
|
| 81 | 166 |
! write(*,*) "Coucou 4-mixed" |
| 82 | 167 |
Call Mixed2Cart(Nat,indzmat,geom,GeomCart) |
| ... | ... | |
| 91 | 176 |
! write(*,*) "Coucou 4-hybrid" |
| 92 | 177 |
GeomCart=reshape(Geom,(/Nat,3/)) |
| 93 | 178 |
CASE DEFAULT |
| 94 |
WRITE (*,*) "Coord=",TRIM(COORD)," not recognized. Stop" |
|
| 95 |
STOP |
|
| 179 |
WRITE (*,*) "Coord=",TRIM(COORD)," not recognized. Stop"
|
|
| 180 |
STOP
|
|
| 96 | 181 |
END SELECT |
| 182 |
!WRITE(*,*) Nat |
|
| 183 |
!WRITE(*,*) 'GeomCart in Egrad_baker' |
|
| 184 |
!DO I=1,Nat |
|
| 185 |
! WRITE(*,'(1X,A,3(1X,F12.8))') AtName(I),GeomCart(I,1:3) |
|
| 186 |
!END DO |
|
| 97 | 187 |
|
| 98 | 188 |
SELECT CASE (Prog) |
| 99 |
CASE ('GAUSSIAN')
|
|
| 189 |
CASE ('GAUSSIAN')
|
|
| 100 | 190 |
! we call the Gaussian routine. |
| 101 | 191 |
! it will return the energy and the gradient in cartesian coordinates. |
| 102 | 192 |
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT) |
| 103 |
Call egrad_gaussian(E,GeomCart,GradCart)
|
|
| 104 |
CASE ('MOPAC')
|
|
| 193 |
Call egrad_gaussian(E,GeomCart,GradCart) |
|
| 194 |
CASE ('MOPAC')
|
|
| 105 | 195 |
! we call the Mopac routine. |
| 106 | 196 |
! it will return the energy and the gradient in cartesian coordinates. |
| 107 | 197 |
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT) |
| 108 |
Call egrad_mopac(E,GeomCart,GradCart)
|
|
| 109 |
CASE ('VASP')
|
|
| 110 |
Call egrad_vasp(E,Geomcart,GradCart)
|
|
| 111 |
CASE ('TURBOMOLE')
|
|
| 112 |
Call egrad_turbomole(E,Geomcart,GradCart)
|
|
| 113 |
CASE ('EXT')
|
|
| 114 |
Call egrad_ext(E,Geomcart,GradCart)
|
|
| 115 |
CASE ('TEST')
|
|
| 116 |
Call egrad_test(Nat,E,Geomcart,GradCart)
|
|
| 117 |
CASE ('CHAMFRE')
|
|
| 118 |
Call egrad_chamfre(Nat,E,Geomcart,GradCart)
|
|
| 198 |
Call egrad_mopac(E,GeomCart,GradCart) |
|
| 199 |
CASE ('VASP')
|
|
| 200 |
Call egrad_vasp(E,Geomcart,GradCart) |
|
| 201 |
CASE ('TURBOMOLE')
|
|
| 202 |
Call egrad_turbomole(E,Geomcart,GradCart) |
|
| 203 |
CASE ('EXT')
|
|
| 204 |
Call egrad_ext(E,Geomcart,GradCart) |
|
| 205 |
CASE ('TEST')
|
|
| 206 |
Call egrad_test(Nat,E,Geomcart,GradCart) |
|
| 207 |
CASE ('CHAMFRE')
|
|
| 208 |
Call egrad_chamfre(Nat,E,Geomcart,GradCart) |
|
| 119 | 209 |
END SELECT |
| 120 | 210 |
if (debug) THEN |
| 121 | 211 |
WRITE(*,*) 'DBG EGRAD, GradCart read' |
| ... | ... | |
| 128 | 218 |
|
| 129 | 219 |
! If (PROG=="VASP") STOP |
| 130 | 220 |
|
| 221 |
|
|
| 131 | 222 |
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid |
| 132 | 223 |
! but we have to convert it into Zmat if COORD=Zmat |
| 133 | 224 |
SELECT CASE (COORD) |
| ... | ... | |
| 162 | 253 |
END DO |
| 163 | 254 |
DEALLOCATE(GradTmp) |
| 164 | 255 |
CASE ('BAKER')
|
| 165 |
WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad_baker.f90, instead of Egrad.f90. Stop" |
|
| 256 |
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid |
|
| 257 |
! but we have to convert it into internal coordinates if COORD=Baker |
|
| 258 |
!!!!!!!!!!!!!!!!!!!! |
|
| 259 |
! |
|
| 260 |
! PFL Debug |
|
| 261 |
! |
|
| 262 |
!!!!!!!!!!!!!!!!!!!! |
|
| 263 |
if (debugPFL) THEN |
|
| 264 |
WRITE(*,*) "DBG PFL Egrad_Baker... computing dzdc" |
|
| 265 |
if (.not.ALLOCATED(indzmat)) THEN |
|
| 266 |
ALLOCATE(indzmat(Nat,5)) |
|
| 267 |
IndZmat=0 |
|
| 268 |
Indzmat(1,1)=1 |
|
| 269 |
IndZmat(2,1)=2 |
|
| 270 |
IndZmat(2,2)=1 |
|
| 271 |
IndZmat(3,1)=3 |
|
| 272 |
IndZmat(3,2)=2 |
|
| 273 |
IndZmat(3,3)=1 |
|
| 274 |
IndZmat(4,1)=4 |
|
| 275 |
IndZmat(4,2)=3 |
|
| 276 |
IndZmat(4,3)=2 |
|
| 277 |
IndZmat(4,4)=1 |
|
| 278 |
END IF |
|
| 279 |
IF (.NOT.ALLOCATED(DzDc)) THEN |
|
| 280 |
ALLOCATE(DzDc(3,nat,3,Nat)) |
|
| 281 |
END IF |
|
| 282 |
if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat)) |
|
| 283 |
DO I=1,Nat |
|
| 284 |
GeomTmp(1,I)=GeomCart(OrderInv(I),1) |
|
| 285 |
GeomTmp(2,I)=GeomCart(OrderInv(i),2) |
|
| 286 |
GeomTmp(3,I)=GeomCart(OrderInv(i),3) |
|
| 287 |
END DO |
|
| 288 |
|
|
| 289 |
DzDc=0.d0 |
|
| 290 |
CALL CalcBMat_int (nat, GeomTmp, indzmat, dzdc, massat,atome) |
|
| 291 |
|
|
| 292 |
END IF ! if (debugPFL) |
|
| 293 |
!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 294 |
! Debugging PFL |
|
| 295 |
!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 296 |
|
|
| 297 |
WRITE(*,*) "BTransInv_local Trans that is used originally" |
|
| 298 |
DO J=1,3*Nat |
|
| 299 |
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) |
|
| 300 |
END DO |
|
| 301 |
|
|
| 302 |
WRITE(*,*) "Calculating actual values using GeomCart" |
|
| 303 |
if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat)) |
|
| 304 |
GeomTmp(1,:)=GeomCart(1:Nat,1) |
|
| 305 |
GeomTmp(2,:)=GeomCart(1:Nat,2) |
|
| 306 |
GeomTmp(3,:)=GeomCart(1:Nat,3) |
|
| 307 |
|
|
| 308 |
! Computing B^prim: |
|
| 309 |
FAloBBT=.NOT.ALLOCATED(BBT) |
|
| 310 |
FAloBBTInv=.NOT.ALLOCATED(BBT_inv) |
|
| 311 |
FAloBPrimT=.NOT.ALLOCATED(BprimT) |
|
| 312 |
if (FAloBPrimT) ALLOCATE(Bprimt(3*Nat,NPrim)) |
|
| 313 |
if (FAloBBT) ALLOCATE(BBt(NCoord,NCoord)) |
|
| 314 |
if (FAloBBTInv) ALLOCATE(BBt_inv(NCoord,NCoord)) |
|
| 315 |
BprimT=0.d0 |
|
| 316 |
ScanCoord=>Coordinate |
|
| 317 |
I=0 |
|
| 318 |
DO WHILE (Associated(ScanCoord%next)) |
|
| 319 |
I=I+1 |
|
| 320 |
SELECT CASE (ScanCoord%Type) |
|
| 321 |
CASE ('BOND')
|
|
| 322 |
CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,GeomTmp,BprimT(1,I)) |
|
| 323 |
CASE ('ANGLE')
|
|
| 324 |
CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,GeomTmp,BprimT(1,I)) |
|
| 325 |
CASE ('DIHEDRAL')
|
|
| 326 |
CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,GeomTmp,BprimT(1,I)) |
|
| 327 |
END SELECT |
|
| 328 |
ScanCoord => ScanCoord%next |
|
| 329 |
END DO |
|
| 330 |
|
|
| 331 |
if (debug) THEN |
|
| 332 |
WRITE(*,*) "Bprim " |
|
| 333 |
DO J=1,3*Nat |
|
| 334 |
WRITE(*,'(50(1X,F12.6))') BprimT(J,:) |
|
| 335 |
END DO |
|
| 336 |
END IF |
|
| 337 |
|
|
| 338 |
BMat_BakerT = 0.d0 |
|
| 339 |
DO I=1,NCoord |
|
| 340 |
DO J=1,NPrim |
|
| 341 |
!BprimT is transpose of B^prim. |
|
| 342 |
!B = UMat^T B^prim, B^T = (B^prim)^T UMat |
|
| 343 |
BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I) |
|
| 344 |
END DO |
|
| 345 |
END DO |
|
| 346 |
|
|
| 347 |
BBT = 0.d0 |
|
| 348 |
DO I=1,NCoord |
|
| 349 |
DO J=1,3*Nat ! BBT(:,I) forms BB^T |
|
| 350 |
BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I) |
|
| 351 |
END DO |
|
| 352 |
END DO |
|
| 353 |
|
|
| 354 |
BBT_inv = 0.d0 |
|
| 355 |
|
|
| 356 |
Call GenInv(NCoord,BBT,BBT_inv,NCoord) |
|
| 357 |
|
|
| 358 |
! ALLOCATE(V(NCoord,NCoord)) |
|
| 359 |
! tol = 1e-10 |
|
| 360 |
! ! BBT is destroyed by GINVSE |
|
| 361 |
! Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT) |
|
| 362 |
! DEALLOCATE(V) |
|
| 363 |
! IF(Fail) Then |
|
| 364 |
! WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90" |
|
| 365 |
! STOP |
|
| 366 |
! END IF |
|
| 367 |
|
|
| 368 |
!Print *, 'BBT_inv=' |
|
| 369 |
DO J=1,NCoord |
|
| 370 |
! WRITE(*,'(50(1X,F10.2))') BBT_inv(:,J) |
|
| 371 |
! Print *, BBT_inv(:,J) |
|
| 372 |
END DO |
|
| 373 |
!Print *, 'End of BBT_inv' |
|
| 374 |
|
|
| 375 |
! Calculation of (B^T)^-1 = (BB^T)^-1B: |
|
| 376 |
BTransInv_local = 0.d0 |
|
| 377 |
DO I=1,3*Nat |
|
| 378 |
DO J=1,NCoord |
|
| 379 |
BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J) |
|
| 380 |
END DO |
|
| 381 |
END DO |
|
| 382 |
|
|
| 383 |
Print *, 'BMat_BakerT=' |
|
| 384 |
DO J=1,3*Nat |
|
| 385 |
WRITE(*,'(50(1X,F12.8))') BMat_BakerT(:,J) |
|
| 386 |
END DO |
|
| 387 |
|
|
| 388 |
WRITE(*,*) "BTransInv_local Trans corresponding to GeomCart" |
|
| 389 |
DO J=1,3*Nat |
|
| 390 |
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) |
|
| 391 |
END DO |
|
| 392 |
|
|
| 393 |
if (FAloBPrimT) DEALLOCATE(Bprimt) |
|
| 394 |
if (FAloBBT) DEALLOCATE(BBt) |
|
| 395 |
if (FAloBBTInv) DEALLOCATE(BBt_inv) |
|
| 396 |
|
|
| 397 |
Grad=0.d0 |
|
| 398 |
DO I=1, 3*Nat |
|
| 399 |
! Here Grad is gradient in Baker coordinates. GradCart(3*Nat) |
|
| 400 |
Grad(:) = Grad(:) + BTransInv_local(:,I)*GradCart(I) !BTransInv_local = BTransInv in Opt_Geom.f90 |
|
| 401 |
END DO |
|
| 402 |
!WRITE(*,*) 'In Egrad.f90, L263, Gradient:' |
|
| 403 |
!WRITE(*,'(12(1X,F8.4))') Grad(1:NCoord) |
|
| 404 |
|
|
| 405 |
|
|
| 166 | 406 |
CASE ('MIXED')
|
| 167 | 407 |
ALLOCATE(GradTmp(3*Nat)) |
| 168 | 408 |
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed" |
| src/Makefile (revision 8) | ||
|---|---|---|
| 5 | 5 |
# You might also have to edit the location of |
| 6 | 6 |
# some libraries (like mkl for ifort) |
| 7 | 7 |
|
| 8 |
Machine=arqg
|
|
| 8 |
Machine=gfortran
|
|
| 9 | 9 |
|
| 10 | 10 |
########################################### |
| 11 | 11 |
# # |
| ... | ... | |
| 148 | 148 |
|
| 149 | 149 |
ifeq ($(Machine),arqg) |
| 150 | 150 |
# Flags for arq GFortran |
| 151 |
COMP=gfortran -fbacktrace -fbounds-check -pedantic -std=gnu
|
|
| 151 |
COMP=gfortran -fbacktrace -fbounds-check -Wall
|
|
| 152 | 152 |
F90=${COMP}
|
| 153 | 153 |
F77=${F90}
|
| 154 | 154 |
#F90=${COMP}
|
| ... | ... | |
| 240 | 240 |
ConvertBakerInternal_cart.f90 \ |
| 241 | 241 |
ConvertZmat_cart_3.f90 \ |
| 242 | 242 |
Egrad.f90 \ |
| 243 |
Egrad_baker.f90 \ |
|
| 244 | 243 |
EgradPath.f90 \ |
| 245 | 244 |
egrad_gaussian.f90 \ |
| 246 | 245 |
egrad_mopac.f90 \ |
| src/Path.f90 (revision 8) | ||
|---|---|---|
| 247 | 247 |
! these are used to read temporary coordinates |
| 248 | 248 |
REAL(KREAL) :: xtmp,ytmp,ztmp |
| 249 | 249 |
|
| 250 |
LOGICAL, ALLOCATABLE :: FTmp(:) ! Nat |
|
| 251 | 250 |
LOGICAL :: FFrozen,FCart |
| 252 | 251 |
|
| 253 | 252 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| src/EgradPath.f90 (revision 8) | ||
|---|---|---|
| 25 | 25 |
REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) |
| 26 | 26 |
REAL(KREAL), ALLOCATABLE :: x_k(:), y_k(:), z_k(:) |
| 27 | 27 |
REAL(KREAL), ALLOCATABLE :: GeomCart(:,:) !(Nat,3) |
| 28 |
REAL(KREAL), ALLOCATABLE :: GeomOld_dummy(:),GeomCart_old_dummy(:,:) |
|
| 29 | 28 |
REAL(KREAL) :: E |
| 30 | 29 |
LOGICAL :: Debug |
| 31 | 30 |
LOGICAL, SAVE :: First=.TRUE. |
| ... | ... | |
| 39 | 38 |
CHARACTER(*), intent(in) :: string |
| 40 | 39 |
logical :: isValid |
| 41 | 40 |
END function VALID |
| 41 |
|
|
| 42 |
subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old) |
|
| 43 |
|
|
| 44 |
! This routines calculates the energy E and the gradient Grad of |
|
| 45 |
! a molecule with Geometry Geom (may be in internal coordinates), |
|
| 46 |
! using for now, either Gaussian or Ext, more general later. |
|
| 47 |
|
|
| 48 |
use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, & |
|
| 49 |
prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, & |
|
| 50 |
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & |
|
| 51 |
, BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & |
|
| 52 |
, Order,OrderInv, XPrimitiveF |
|
| 53 |
|
|
| 54 |
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord) |
|
| 55 |
! allocated in Path.f90 |
|
| 56 |
|
|
| 57 |
use Io_module |
|
| 58 |
|
|
| 59 |
! Energy (calculated if F300K=.F., else estimated) |
|
| 60 |
REAL(KREAL), INTENT (OUT) :: E |
|
| 61 |
! NCoord: Number of the degrees of freedom |
|
| 62 |
! IGeom: index of the geometry. |
|
| 63 |
INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt |
|
| 64 |
! Geometry at which gradient is calculated (cf Factual also): |
|
| 65 |
REAL(KREAL), INTENT (INOUT) :: Geom(NCoord) |
|
| 66 |
! Gradient calculated at Geom geometry: |
|
| 67 |
REAL(KREAL), INTENT (OUT) :: Grad(NCoord) |
|
| 68 |
! Cartesian geometry corresponding to (Internal Geometry) Geom: |
|
| 69 |
REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3) |
|
| 70 |
!!! Optional, just for geometry optimization with Baker coordinates |
|
| 71 |
REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3) |
|
| 72 |
REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord) |
|
| 73 |
! FOptGeom is a flag indicating if we are doing a geom optimization |
|
| 74 |
! it can be omitted so that we use a local flag: Flag_Opt_Geom |
|
| 75 |
LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom |
|
| 76 |
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F |
|
| 77 |
LOGICAL :: Flag_Opt_Geom |
|
| 78 |
|
|
| 79 |
END subroutine Egrad |
|
| 80 |
|
|
| 42 | 81 |
END INTERFACE |
| 43 | 82 |
|
| 44 | 83 |
debug=valid('EGradPath')
|
| ... | ... | |
| 48 | 87 |
ALLOCATE(GeomTmp(NCoord)) |
| 49 | 88 |
ALLOCATE(x(Nat),y(Nat),z(Nat)) |
| 50 | 89 |
ALLOCATE(x_k(Nat),y_k(Nat),z_k(Nat)) |
| 51 |
ALLOCATE(GeomOld_dummy(NCoord)) |
|
| 52 |
ALLOCATE(GeomCart_old_dummy(Nat,3)) |
|
| 53 | 90 |
|
| 54 | 91 |
|
| 55 | 92 |
IF (RunMode=="PARA") THEN ! matches at L315 |
| ... | ... | |
| 608 | 645 |
ELSE ! matches IF (RunMode=="PARA") THEN |
| 609 | 646 |
! We will launch all calculations sequentially |
| 610 | 647 |
ALLOCATE(GradTmp(NCoord)) |
| 611 |
GeomOld_dummy=0.d0 ! Internal coordinates |
|
| 612 |
GeomCart_old_dummy=0.d0 |
|
| 613 | 648 |
! We have the new path, we calculate its energy and gradient |
| 614 | 649 |
IGeom0=2 |
| 615 | 650 |
IGeomF=NGeomF-1 |
| ... | ... | |
| 664 | 699 |
WRITE(*,'(12(1X,F6.3))') GeomCart |
| 665 | 700 |
END IF |
| 666 | 701 |
|
| 667 |
IF (COORD.EQ.'BAKER') THEN |
|
| 668 | 702 |
! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT). |
| 669 | 703 |
! GeomCart has INTENT(OUT) |
| 670 |
Call EGrad_baker(E,GeomTmp,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, & |
|
| 671 |
GeomOld_dummy,GeomCart_old_dummy) |
|
| 672 |
ELSE |
|
| 673 |
Call EGrad(E,GeomTmp,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom) !GeomTmp=IntCoordF(IGeom,:) |
|
| 674 |
END IF |
|
| 704 |
Call EGrad(E,GeomTmp,GradTmp,NCoord,IGeom,IOpt,GeomCart) !GeomTmp=IntCoordF(IGeom,:) |
|
| 675 | 705 |
|
| 676 | 706 |
! Egrad calls ConvertBakerInternal_cart to convert GeomTmp=IntCoordF(IGeom,:) into |
| 677 | 707 |
! GeomCart (Cartesian Coordinates) so that energy and gradients can be calculated |
| ... | ... | |
| 694 | 724 |
DEALLOCATE(GeomTmp) |
| 695 | 725 |
DEALLOCATE(x,y,z) |
| 696 | 726 |
DEALLOCATE(x_k,y_k,z_k) |
| 697 |
DEALLOCATE(GeomOld_dummy,GeomCart_old_dummy) |
|
| 698 | 727 |
if (debug) Call header('Exiting EgradPath')
|
| 699 | 728 |
RETURN |
| 700 | 729 |
999 WRITE(*,*) "EgradPath : We should not be here !" |
| src/Opt_Geom.f90 (revision 8) | ||
|---|---|---|
| 18 | 18 |
LOGICAL, INTENT(IN) :: Flag_Opt_Geom |
| 19 | 19 |
|
| 20 | 20 |
! As this subroutine is here for debugging, debug=T ! |
| 21 |
LOGICAL, PARAMETER :: debug=.TRUE.
|
|
| 21 |
LOGICAL :: debug
|
|
| 22 | 22 |
LOGICAL :: Fini |
| 23 | 23 |
LOGICAL, SAVE :: FRST=.TRUE. |
| 24 | 24 |
|
| ... | ... | |
| 38 | 38 |
REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis |
| 39 | 39 |
REAL(KREAL) :: NormStep, FactStep, HP |
| 40 | 40 |
|
| 41 |
INTERFACE |
|
| 42 |
function valid(string) result (isValid) |
|
| 43 |
CHARACTER(*), intent(in) :: string |
|
| 44 |
logical :: isValid |
|
| 45 |
END function VALID |
|
| 46 |
|
|
| 47 |
subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old) |
|
| 48 |
|
|
| 49 |
! This routines calculates the energy E and the gradient Grad of |
|
| 50 |
! a molecule with Geometry Geom (may be in internal coordinates), |
|
| 51 |
! using for now, either Gaussian or Ext, more general later. |
|
| 52 |
|
|
| 53 |
use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, & |
|
| 54 |
prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, & |
|
| 55 |
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & |
|
| 56 |
, BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & |
|
| 57 |
, Order,OrderInv, XPrimitiveF |
|
| 58 |
|
|
| 59 |
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord) |
|
| 60 |
! allocated in Path.f90 |
|
| 61 |
|
|
| 62 |
use Io_module |
|
| 63 |
|
|
| 64 |
! Energy (calculated if F300K=.F., else estimated) |
|
| 65 |
REAL(KREAL), INTENT (OUT) :: E |
|
| 66 |
! NCoord: Number of the degrees of freedom |
|
| 67 |
! IGeom: index of the geometry. |
|
| 68 |
INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt |
|
| 69 |
! Geometry at which gradient is calculated (cf Factual also): |
|
| 70 |
REAL(KREAL), INTENT (INOUT) :: Geom(NCoord) |
|
| 71 |
! Gradient calculated at Geom geometry: |
|
| 72 |
REAL(KREAL), INTENT (OUT) :: Grad(NCoord) |
|
| 73 |
! Cartesian geometry corresponding to (Internal Geometry) Geom: |
|
| 74 |
REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3) |
|
| 75 |
!!! Optional, just for geometry optimization with Baker coordinates |
|
| 76 |
REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3) |
|
| 77 |
REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord) |
|
| 78 |
! FOptGeom is a flag indicating if we are doing a geom optimization |
|
| 79 |
! it can be omitted so that we use a local flag: Flag_Opt_Geom |
|
| 80 |
LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom |
|
| 81 |
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F |
|
| 82 |
LOGICAL :: Flag_Opt_Geom |
|
| 83 |
|
|
| 84 |
END subroutine Egrad |
|
| 85 |
|
|
| 86 |
END INTERFACE |
|
| 87 |
|
|
| 88 |
debug=valid('OptGeom')
|
|
| 89 |
|
|
| 41 | 90 |
E=0. |
| 42 | 91 |
|
| 43 |
if (debug) WRITE(*,*) "=================== Geom Optimization =================="
|
|
| 92 |
if (debug) Call Header("Entering Geom Optimization")
|
|
| 44 | 93 |
|
| 45 | 94 |
ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord)) |
| 46 | 95 |
ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord)) |
| ... | ... | |
| 190 | 239 |
! GeomCart has INTENT(OUT) |
| 191 | 240 |
! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable |
| 192 | 241 |
GeomRef=GeomOld |
| 193 |
Call EGrad_baker(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, &
|
|
| 242 |
Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, & |
|
| 194 | 243 |
GeomRef,GeomCart_old) |
| 195 | 244 |
GeomCart_old=GeomCart |
| 196 | 245 |
ELSE |
| ... | ... | |
| 378 | 427 |
END IF |
| 379 | 428 |
|
| 380 | 429 |
DEALLOCATE(GeomCart) |
| 381 |
if (debug) WRITE(*,*) "================== Geom Optimization Over=================" |
|
| 430 |
|
|
| 431 |
if (debug) Call Header("Geom Optimization Over")
|
|
| 432 |
|
|
| 382 | 433 |
END SUBROUTINE Opt_geom |
Formats disponibles : Unified diff