Statistics
| Revision:

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 re-orient 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 4-mixed"
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*Iat-2:3*Iat)
238
        END DO
239
     END IF
240
  CASE ('HYBRID')
241
!  write(*,*) "Coucou 4-hybrid"
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*I-2: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*i-2)
323
        Grad(Idx+1)=GradTmp(3*i-1)
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 = 1e-10
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*I-2: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*i-2)
519
        Grad(Idx+1)=GradTmp(3*i-1)
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