Révision 8 src/Egrad.f90
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 4mixed" 
82  167 
Call Mixed2Cart(Nat,indzmat,geom,GeomCart) 
...  ...  
91  176 
! write(*,*) "Coucou 4hybrid" 
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 = 1e10 

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" 
Formats disponibles : Unified diff