root / src / Egrad.f90 @ 11
History  View  Annotate  Download (17.3 kB)
1 
subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,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,Atome,massat, & 
8 
prog,NCart,XyzGeomF, XyzGeomI, renum, & 
9 
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & 
10 
, BprimT,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & 
11 
,Order, OrderInv, XPrimitiveF,FPBC,XGeomRefPBC,YGeomRefPBC,ZGeomRefPBC 
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 
! Gradient calculated at Geom geometry: 
27 
REAL(KREAL), INTENT (OUT) :: Grad(NCoord) 
28 
! Cartesian geometry corresponding to (Internal Geometry) Geom: 
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 
38  
39 
! ====================================================================== 
40  
41 
LOGICAL :: debug 
42 
REAL(KREAL), ALLOCATABLE :: GradTmp(:) 
43 
REAL(KREAL), ALLOCATABLE :: GradCart(:) 
44 
REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) 
45 
REAL(KREAL) :: Pi 
46 
REAL(KREAL) :: MRot(3,3), Rmsd 
47  
48 
INTEGER(KINT) :: iat, i, j, IBeg,Idx 
49  
50 
REAL(KREAL), ALLOCATABLE :: XPrimRef(:), XPrim(:) ! NPrim 
51  
52 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
53 
! 
54 
!!!!!!!! PFL Debug 
55 
! 
56 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
57 
REAL(KREAL), ALLOCATABLE :: GeomTmp(:,:) !3,Nat 
58 
LOGICAL :: FALOBBT,FALOBPrimt,FAloBBTInv 
59 
LOGICAL :: DebugPFL 
60  
61 
! ====================================================================== 
62 
INTERFACE 
63 
function valid(string) result (isValid) 
64 
CHARACTER(*), intent(in) :: string 
65 
logical :: isValid 
66 
END function VALID 
67 
END INTERFACE 
68 
! ====================================================================== 
69  
70 
Pi=dacos(1.0d0) 
71  
72 
if (present(FOptGeom)) THEN 
73 
Flag_Opt_Geom=FOptGeom 
74 
ELSE 
75 
Flag_Opt_Geom=.FALSE. 
76 
END IF 
77  
78 
debug=valid('EGRAD') 
79 
debugPFL=valid('bakerPFL') 
80 
if (debug) Call Header("Entering Egrad") 
81  
82 
if (debug) Print *, 'EGrad.f90, L73, IOpt=', IOpt, ' Flag_Opt_Geom=',Flag_Opt_Geom 
83  
84 
Grad=0. 
85 
E=0. 
86  
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 
! In order to avoid problem with numbering and linear angles/diehedral, we convert the 
95 
! zmat into cartesian coordinates 
96 
Call Int2Cart(nat,indzmat,Geom,GeomCart) 
97 
! In case of PBC calculations, we have to reorient the geometry onto the user initial geometry 
98 
If (FPBC) THEN 
99 
! we align this geometry on the initial one 
100 
if (debug) THEN 
101 
WRITe(*,*) "We are orientating..." 
102 
WRITE(*,*) "Geom before orientation" 
103 
WRITE(*,*) Nat 
104 
WRITE(*,*) "" 
105 
DO I=1,Nat 
106 
If (renum) Iat=Order(I) 
107 
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomCart(I,1),GeomCart(I,2),GeomCart(I,3) 
108 
END DO 
109 
END IF 
110  
111 
IF (Renum) THEN 
112 
DO I=1,Nat 
113 
Iat=Order(I) 
114 
X(I)=GeomCart(Iat,1) 
115 
Y(I)=GeomCart(Iat,2) 
116 
Z(I)=GeomCart(Iat,3) 
117 
END DO 
118 
ELSE 
119 
X(1:Nat)=GeomCart(1:Nat,1) 
120 
Y(1:Nat)=GeomCart(1:Nat,2) 
121 
Z(1:Nat)=GeomCart(1:Nat,3) 
122 
END IF 
123  
124 
if (debug) then 
125 
WRITE(*,*) "Geom before orientation after reorderring" 
126 
WRITE(*,*) Nat 
127 
WRITE(*,*) "" 
128 
DO I=1,Nat 
129 
! Iat=I 
130 
! If (Renum) Iat=Order(I) 
131 
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),X(I),Y(I),Z(I) 
132 
END DO 
133 
WRITE(*,*) "Ref Geom" 
134 
WRITE(*,*) Nat 
135 
WRITE(*,*) "" 
136 
DO I=1,Nat 
137 
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),XGeomRefPBC(I),YGeomRefPBC(I),ZGeomRefPBC(I) 
138 
END DO 
139 
END IF 
140  
141 
Call CalcRmsd(Nat,XGeomRefPBC,YGeomRefPBC,ZGeomRefPBC, & 
142 
X,Y,Z, MRot, Rmsd,.TRUE.,.TRUE.) 
143  
144 
If (DEBUG) THEN 
145 
WRITE(*,*) "Geom AFTER orientation" 
146 
WRITE(*,*) Nat 
147 
WRITE(*,*) "" 
148 
DO I=1,Nat 
149 
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),X(I),Y(I),Z(I) 
150 
END DO 
151 
END IF 
152 

153 
If (Renum) THEN 
154 
Do I=1,Nat 
155 
Iat=orderInv(I) 
156 
GeomCart(Iat,1)=X(I) 
157 
GeomCart(Iat,2)=Y(I) 
158 
GeomCart(Iat,3)=Z(I) 
159 
END DO 
160 
END IF 
161 
END IF 
162  
163 
CASE ('BAKER') 
164 
XPrimRef=XPrimitiveF(IGeom,:) 
165 
IF (Flag_Opt_Geom) THEN 
166 
IF (IOpt .EQ. 1) THEN 
167 
GeomCart(:,1) = XyzGeomI(IGeom,1,:) 
168 
GeomCart(:,2) = XyzGeomI(IGeom,2,:) 
169 
GeomCart(:,3) = XyzGeomI(IGeom,3,:) 
170 
ELSE 
171 
if (debug) WRITE(*,*) 'Before Call ConvertBakerInternal_cart; L125, Egrad.f90' 
172 
if (.NOT.present(GeomOld)) THEN 
173 
WRITE(*,*) "ERROR: GeomOld not passed as an argument" 
174 
STOP 
175 
END IF 
176 
WRITE(*,*) 'GeomOld=' 
177 
WRITE(*,'(12(1X,F6.3))') GeomOld 
178 
WRITE(*,*) 'Geom=' 
179 
WRITE(*,'(12(1X,F6.3))') Geom 
180 
if (.NOT.present(GeomCart_Old)) THEN 
181 
WRITE(*,*) "ERROR: GeomCart_Old not passed as an argument" 
182 
STOP 
183 
END IF 
184  
185 
WRITE(*,*) 'GeomCart_old=' 
186 
WRITE(*,'(20(1X,F6.3))') GeomCart_old(1:Nat,1:3) 
187 
Call ConvertBakerInternal_cart(GeomOld,Geom,GeomCart_old(:,1), & 
188 
GeomCart_old(:,2),GeomCart_old(:,3),x,y,z,XPrim,XPrimRef) 
189 
GeomCart(:,1) = x(:) 
190 
GeomCart(:,2) = y(:) 
191 
GeomCart(:,3) = z(:) 
192 
END IF 
193 
ELSE 
194 
IF (IOpt .EQ. 0) THEN 
195 
GeomCart(:,1) = XyzGeomF(IGeom,1,:) 
196 
GeomCart(:,2) = XyzGeomF(IGeom,2,:) 
197 
GeomCart(:,3) = XyzGeomF(IGeom,3,:) 
198 
ELSE 
199 
! IntCoordF(NGeomF,NCoord), XyzGeomF(NGeomF,3,Nat) 
200 
! Geom has to be converted into x,y,z 
201 
! GeomOld_all and XyzGeomF are old geometry in internal and cartesian coords resp. 
202 
DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart() 
203 
BTransInv_local(I,:) = BTransInvF(IGeom,I,:) 
204 
UMat_local(:,I) = UMatF(IGeom,:,I) 
205 
END DO 
206 
if (debug) Print *, 'I am before Call ConvertBakerInternal_cart; L160, Egrad.f90' 
207 
WRITE(*,*) 'GeomOld=' 
208 
WRITE(*,'(12(1X,F6.3))') GeomOld_all(IGeom,:) 
209 
WRITE(*,*) 'Geom=' 
210 
WRITE(*,'(12(1X,F6.3))') Geom 
211 
WRITE(*,*) 'DBG Egrad L165 GeomCart_old=' 
212 
DO I=1,Nat 
213 
WRITE(*,'(3(1X,F12.8))') XyzGeomF(IGeom,:,I) 
214 
END DO 
215  
216 
Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),Geom,XyzGeomF(IGeom,1,:), & 
217 
XyzGeomF(IGeom,2,:),XyzGeomF(IGeom,3,:),x,y,z,XPrim,XPrimRef) 
218 
if (debug) Print *, 'I am after Call ConvertBakerInternal_cart; L111, Egrad.f90' 
219 
DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart() 
220 
BTransInvF(IGeom,I,:) = BTransInv_local(I,:) 
221 
END DO 
222 
! GeomCart is actually XyzCoord and has INTENT(OUT) in this subroutine. 
223 
GeomCart(:,1) = x(:) 
224 
GeomCart(:,2) = y(:) 
225 
GeomCart(:,3) = z(:) 
226 
END IF ! matches IF (IOpt .EQ. 0) THEN 
227 
END IF 
228 
CASE ('MIXED') 
229 
! write(*,*) "Coucou 4mixed" 
230 
Call Mixed2Cart(Nat,indzmat,geom,GeomCart) 
231 
CASE ('CART') 
232  
233 
GeomCart=reshape(Geom,(/Nat,3/)) 
234 
if (debug) THEN 
235 
WRITE(*,*) "Coord=CART... in egrad" 
236 
DO Iat=1,Nat 
237 
WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),Geom(3*Iat2:3*Iat) 
238 
END DO 
239 
END IF 
240 
CASE ('HYBRID') 
241 
! write(*,*) "Coucou 4hybrid" 
242 
GeomCart=reshape(Geom,(/Nat,3/)) 
243 
CASE DEFAULT 
244 
WRITE (*,*) "Coord=",TRIM(COORD)," not recognized. Stop" 
245 
STOP 
246 
END SELECT 
247 
!WRITE(*,*) Nat 
248 
!WRITE(*,*) 'GeomCart in Egrad_baker' 
249 
!DO I=1,Nat 
250 
! WRITE(*,'(1X,A,3(1X,F12.8))') AtName(I),GeomCart(I,1:3) 
251 
!END DO 
252  
253  
254  
255 
SELECT CASE (Prog) 
256 
CASE ('GAUSSIAN') 
257 
! we call the Gaussian routine. 
258 
! it will return the energy and the gradient in cartesian coordinates. 
259 
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT) 
260 
Call egrad_gaussian(E,GeomCart,GradCart) 
261 
CASE ('MOPAC') 
262 
! we call the Mopac routine. 
263 
! it will return the energy and the gradient in cartesian coordinates. 
264 
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT) 
265 
Call egrad_mopac(E,GeomCart,GradCart) 
266 
CASE ('VASP') 
267 
Call egrad_vasp(E,Geomcart,GradCart) 
268 
CASE ('SIESTA') 
269 
Call egrad_siesta(E,Geomcart,GradCart) 
270 
CASE ('TURBOMOLE') 
271 
Call egrad_turbomole(E,Geomcart,GradCart) 
272 
CASE ('EXT') 
273 
Call egrad_ext(E,Geomcart,GradCart) 
274 
CASE ('TEST') 
275 
Call egrad_test(Nat,E,Geomcart,GradCart) 
276 
CASE ('TEST2D') 
277 
Call egrad_test_2D(Nat,E,Geomcart,GradCart) 
278 
CASE ('CHAMFRE') 
279 
Call egrad_chamfre(Nat,E,Geomcart,GradCart) 
280 
CASE ('LEPS') 
281 
Call egrad_LEPS(Nat,E,Geomcart,GradCart) 
282 
END SELECT 
283 
if (debug) THEN 
284 
WRITE(*,*) 'DBG EGRAD, GradCart read' 
285 
DO I=1,Nat 
286 
Iat=I 
287 
if (renum) Iat=orderInv(I) 
288 
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GradCart(3*I2:3*I) 
289 
END DO 
290 
END IF 
291  
292 
! If (PROG=="VASP") STOP 
293  
294  
295 
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid 
296 
! but we have to convert it into Zmat if COORD=Zmat 
297 
SELECT CASE (COORD) 
298 
CASE ("ZMAT") 
299 
ALLOCATE(GradTmp(3*Nat)) 
300 
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int" 
301 
CALL CalcBMat_int (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome) 
302  
303 
if (debug) WRITE(*,*) "DBG EGRAD Calling ConvGrad_Cart2Int" 
304 
CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp) 
305  
306 
if (debug) WRITE(*,*) "DBG EGRAD Storing Grad" 
307 
Grad(1)=GradTmp(4) 
308 
Grad(2)=GradTmp(7) 
309 
! We might have troubles whith angles that are not in the [0:pi] range because, 
310 
! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian... 
311 
! so that the derivative is not good, and a multiplicative factor should be added... 
312 
! 
313 
! This in fact should be taken care of in B mat calculation... 
314 
! 
315 
IF ((Geom(3).LE.0).OR.(Geom(3).GE.Pi)) THEN 
316 
Grad(3)=1.0d0*GradTmp(8) 
317 
ELSE 
318 
Grad(3)=GradTmp(8) 
319 
END IF 
320 
Idx=4 
321 
DO I=4,Nat 
322 
Grad(Idx)=GradTmp(3*i2) 
323 
Grad(Idx+1)=GradTmp(3*i1) 
324 
Grad(Idx+2)=GradTmp(3*i) 
325 
Idx=Idx+3 
326 
END DO 
327 
DEALLOCATE(GradTmp) 
328 
CASE ('BAKER') 
329 
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid 
330 
! but we have to convert it into internal coordinates if COORD=Baker 
331 
!!!!!!!!!!!!!!!!!!!! 
332 
! 
333 
! PFL Debug 
334 
! 
335 
!!!!!!!!!!!!!!!!!!!! 
336 
if (debugPFL) THEN 
337 
WRITE(*,*) "DBG PFL Egrad_Baker... computing dzdc" 
338 
if (.not.ALLOCATED(indzmat)) THEN 
339 
ALLOCATE(indzmat(Nat,5)) 
340 
IndZmat=0 
341 
Indzmat(1,1)=1 
342 
IndZmat(2,1)=2 
343 
IndZmat(2,2)=1 
344 
IndZmat(3,1)=3 
345 
IndZmat(3,2)=2 
346 
IndZmat(3,3)=1 
347 
IndZmat(4,1)=4 
348 
IndZmat(4,2)=3 
349 
IndZmat(4,3)=2 
350 
IndZmat(4,4)=1 
351 
END IF 
352 
IF (.NOT.ALLOCATED(DzDc)) THEN 
353 
ALLOCATE(DzDc(3,nat,3,Nat)) 
354 
END IF 
355 
if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat)) 
356 
DO I=1,Nat 
357 
GeomTmp(1,I)=GeomCart(OrderInv(I),1) 
358 
GeomTmp(2,I)=GeomCart(OrderInv(i),2) 
359 
GeomTmp(3,I)=GeomCart(OrderInv(i),3) 
360 
END DO 
361  
362 
DzDc=0.d0 
363 
CALL CalcBMat_int (nat, GeomTmp, indzmat, dzdc, massat,atome) 
364  
365 
END IF ! if (debugPFL) 
366 
!!!!!!!!!!!!!!!!!!!!!!!!! 
367 
! Debugging PFL 
368 
!!!!!!!!!!!!!!!!!!!!!!!! 
369  
370 
WRITE(*,*) "BTransInv_local Trans that is used originally" 
371 
DO J=1,3*Nat 
372 
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) 
373 
END DO 
374  
375 
WRITE(*,*) "Calculating actual values using GeomCart" 
376 
if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat)) 
377 
GeomTmp(1,:)=GeomCart(1:Nat,1) 
378 
GeomTmp(2,:)=GeomCart(1:Nat,2) 
379 
GeomTmp(3,:)=GeomCart(1:Nat,3) 
380  
381 
! Computing B^prim: 
382 
FAloBBT=.NOT.ALLOCATED(BBT) 
383 
FAloBBTInv=.NOT.ALLOCATED(BBT_inv) 
384 
FAloBPrimT=.NOT.ALLOCATED(BprimT) 
385 
if (FAloBPrimT) ALLOCATE(Bprimt(3*Nat,NPrim)) 
386 
if (FAloBBT) ALLOCATE(BBt(NCoord,NCoord)) 
387 
if (FAloBBTInv) ALLOCATE(BBt_inv(NCoord,NCoord)) 
388 
BprimT=0.d0 
389 
ScanCoord=>Coordinate 
390 
I=0 
391 
DO WHILE (Associated(ScanCoord%next)) 
392 
I=I+1 
393 
SELECT CASE (ScanCoord%Type) 
394 
CASE ('BOND') 
395 
CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,GeomTmp,BprimT(1,I)) 
396 
CASE ('ANGLE') 
397 
CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,GeomTmp,BprimT(1,I)) 
398 
CASE ('DIHEDRAL') 
399 
CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,GeomTmp,BprimT(1,I)) 
400 
END SELECT 
401 
ScanCoord => ScanCoord%next 
402 
END DO 
403  
404 
if (debug) THEN 
405 
WRITE(*,*) "Bprim " 
406 
DO J=1,3*Nat 
407 
WRITE(*,'(50(1X,F12.6))') BprimT(J,:) 
408 
END DO 
409 
END IF 
410  
411 
BMat_BakerT = 0.d0 
412 
DO I=1,NCoord 
413 
DO J=1,NPrim 
414 
!BprimT is transpose of B^prim. 
415 
!B = UMat^T B^prim, B^T = (B^prim)^T UMat 
416 
BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I) 
417 
END DO 
418 
END DO 
419  
420 
BBT = 0.d0 
421 
DO I=1,NCoord 
422 
DO J=1,3*Nat ! BBT(:,I) forms BB^T 
423 
BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I) 
424 
END DO 
425 
END DO 
426  
427 
BBT_inv = 0.d0 
428  
429 
Call GenInv(NCoord,BBT,BBT_inv,NCoord) 
430  
431 
! ALLOCATE(V(NCoord,NCoord)) 
432 
! tol = 1e10 
433 
! ! BBT is destroyed by GINVSE 
434 
! Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT) 
435 
! DEALLOCATE(V) 
436 
! IF(Fail) Then 
437 
! WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90" 
438 
! STOP 
439 
! END IF 
440  
441 
!Print *, 'BBT_inv=' 
442 
DO J=1,NCoord 
443 
! WRITE(*,'(50(1X,F10.2))') BBT_inv(:,J) 
444 
! Print *, BBT_inv(:,J) 
445 
END DO 
446 
!Print *, 'End of BBT_inv' 
447  
448 
! Calculation of (B^T)^1 = (BB^T)^1B: 
449 
BTransInv_local = 0.d0 
450 
DO I=1,3*Nat 
451 
DO J=1,NCoord 
452 
BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J) 
453 
END DO 
454 
END DO 
455  
456 
Print *, 'BMat_BakerT=' 
457 
DO J=1,3*Nat 
458 
WRITE(*,'(50(1X,F12.8))') BMat_BakerT(:,J) 
459 
END DO 
460  
461 
WRITE(*,*) "BTransInv_local Trans corresponding to GeomCart" 
462 
DO J=1,3*Nat 
463 
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) 
464 
END DO 
465  
466 
if (FAloBPrimT) DEALLOCATE(Bprimt) 
467 
if (FAloBBT) DEALLOCATE(BBt) 
468 
if (FAloBBTInv) DEALLOCATE(BBt_inv) 
469  
470 
Grad=0.d0 
471 
DO I=1, 3*Nat 
472 
! Here Grad is gradient in Baker coordinates. GradCart(3*Nat) 
473 
Grad(:) = Grad(:) + BTransInv_local(:,I)*GradCart(I) !BTransInv_local = BTransInv in Opt_Geom.f90 
474 
END DO 
475 
!WRITE(*,*) 'In Egrad.f90, L263, Gradient:' 
476 
!WRITE(*,'(12(1X,F8.4))') Grad(1:NCoord) 
477  
478  
479 
CASE ('MIXED') 
480 
ALLOCATE(GradTmp(3*Nat)) 
481 
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed" 
482 
CALL CalcBMat_mixed (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome) 
483  
484 
if (Debug) THEN 
485 
WRITE(*,*) "DzDc" 
486 
DO I=1,Nat 
487 
DO J=1,3 
488 
WRITE(*,'(50(1X,F12.6))') dzdc(J,I,:,:) 
489 
END DO 
490 
END DO 
491 
END IF 
492 

493  
494 
CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp) 
495 
DO I=1,Nat 
496 
WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I2:3*I) 
497 
END DO 
498 
SELECT CASE (NCART) 
499 
CASE (1) 
500 
Grad(1:3)=GradTmp(1:3) 
501 
Grad(4)=GradTmp(4) 
502 
Grad(5)=GradTmp(7) 
503 
Grad(6)=GradTmp(8) 
504 
Idx=7 
505 
IBeg=4 
506 
CASE(2) 
507 
Grad(1:3)=GradTmp(1:3) 
508 
Grad(4:6)=GradTmp(4:6) 
509 
Grad(7)=GradTmp(7) 
510 
Grad(8)=GradTmp(8) 
511 
Idx=9 
512 
IBeg=4 
513 
CASE DEFAULT 
514 
Idx=1 
515 
IBeg=1 
516 
END SELECT 
517 
DO I=IBeg,Nat 
518 
Grad(Idx)=GradTmp(3*i2) 
519 
Grad(Idx+1)=GradTmp(3*i1) 
520 
Grad(Idx+2)=GradTmp(3*i) 
521 
Idx=Idx+3 
522 
END DO 
523 
DEALLOCATE(GradTmp) 
524 
CASE ("CART","HYBRID") 
525 
Grad=GradCart 
526 
CASE DEFAULT 
527 
WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in Egradients. STOP" 
528 
STOP 
529 
END SELECT 
530  
531 
DEALLOCATE(GradCart) 
532 
DEALLOCATE(x,y,z) 
533  
534 
IF (Debug) WRITE(*,*) 'DBG EGRAD grad',grad(1:NCoord) 
535 
if (debug) Call Header("Egrad Over") 
536  
537 
!999 CONTINUE 
538 
!if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!' 
539 
!STOP 
540 
! ====================================================================== 
541 
end subroutine Egrad 
542 