Statistiques
| Révision :

root / src / Egrad.f90 @ 8

Historique | Voir | Annoter | Télécharger (15,52 ko)

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
      , 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
  ! 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

    
47
  INTEGER(KINT) :: iat, i, j, IBeg,Idx
48

    
49
  REAL(KREAL), ALLOCATABLE :: XPrimRef(:), XPrim(:) ! NPrim
50

    
51
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
52
!
53
!!!!!!!! PFL Debug
54
!
55
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
56
  REAL(KREAL), ALLOCATABLE :: GeomTmp(:,:) !3,Nat
57
  LOGICAL ::  FALOBBT,FALOBPrimt,FAloBBTInv
58
  LOGICAL :: DebugPFL
59

    
60
  ! ======================================================================
61
  INTERFACE
62
     function valid(string) result (isValid)
63
       CHARACTER(*), intent(in) :: string
64
       logical                  :: isValid
65
     END function VALID
66
  END INTERFACE
67
  ! ======================================================================
68

    
69
  Pi=dacos(-1.0d0)
70

    
71
  if (present(FOptGeom)) THEN
72
     Flag_Opt_Geom=FOptGeom
73
  ELSE
74
     Flag_Opt_Geom=.FALSE.
75
  END IF
76

    
77
  debug=valid('EGRAD')
78
  debugPFL=valid('bakerPFL')
79
  if (debug) Call Header("Entering Egrad")
80

    
81
  if (debug) Print *, 'EGrad.f90, L73, IOpt=', IOpt, ' Flag_Opt_Geom=',Flag_Opt_Geom
82

    
83
  Grad=0.
84
  E=0.
85

    
86

    
87
  ALLOCATE(GradCart(3*Nat))
88
  ALLOCATE(x(Nat),y(Nat),z(Nat))
89
  ALLOCATE(XPrimRef(NPrim),XPrim(NPrim))
90

    
91
  SELECT CASE (COORD)
92
  CASE ('ZMAT')
93
     ! In order to avoid problem with numbering and linear angles/diehedral, we convert the
94
     ! zmat into cartesian coordinates
95
     ! A remplacer par Int2Cart :-)
96
     Call Int2Cart(nat,indzmat,Geom,GeomCart)
97
  CASE ('BAKER')
98
     XPrimRef=XPrimitiveF(IGeom,:)
99
     IF (Flag_Opt_Geom) THEN
100
        IF (IOpt .EQ. 1) THEN
101
           GeomCart(:,1) = XyzGeomI(IGeom,1,:)
102
           GeomCart(:,2) = XyzGeomI(IGeom,2,:)
103
           GeomCart(:,3) = XyzGeomI(IGeom,3,:)
104
        ELSE
105
           if (debug) WRITE(*,*) 'Before Call ConvertBakerInternal_cart; L125, Egrad.f90'
106
           if (.NOT.present(GeomOld)) THEN
107
              WRITE(*,*) "ERROR: GeomOld not passed as an argument"
108
              STOP
109
           END IF
110
           WRITE(*,*) 'GeomOld='
111
           WRITE(*,'(12(1X,F6.3))') GeomOld
112
           WRITE(*,*) 'Geom='
113
           WRITE(*,'(12(1X,F6.3))') Geom
114
           if (.NOT.present(GeomCart_Old)) THEN
115
              WRITE(*,*) "ERROR: GeomCart_Old not passed as an argument"
116
              STOP
117
           END IF
118

    
119
           WRITE(*,*) 'GeomCart_old='
120
           WRITE(*,'(20(1X,F6.3))') GeomCart_old(1:Nat,1:3)
121
           Call ConvertBakerInternal_cart(GeomOld,Geom,GeomCart_old(:,1), &
122
                GeomCart_old(:,2),GeomCart_old(:,3),x,y,z,XPrim,XPrimRef)
123
           GeomCart(:,1) = x(:)
124
           GeomCart(:,2) = y(:)
125
           GeomCart(:,3) = z(:)
126
        END IF
127
     ELSE
128
        IF (IOpt .EQ. 0) THEN  
129
           GeomCart(:,1) = XyzGeomF(IGeom,1,:)
130
           GeomCart(:,2) = XyzGeomF(IGeom,2,:)
131
           GeomCart(:,3) = XyzGeomF(IGeom,3,:)
132
        ELSE
133
           ! IntCoordF(NGeomF,NCoord), XyzGeomF(NGeomF,3,Nat)
134
           ! Geom has to be converted into x,y,z
135
           ! GeomOld_all and XyzGeomF are old geometry in internal and cartesian coords resp.
136
           DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart()
137
              BTransInv_local(I,:) = BTransInvF(IGeom,I,:)
138
              UMat_local(:,I) = UMatF(IGeom,:,I)
139
           END DO
140
           if (debug) Print *, 'I am before Call ConvertBakerInternal_cart; L160, Egrad.f90'
141
           WRITE(*,*) 'GeomOld='
142
           WRITE(*,'(12(1X,F6.3))') GeomOld_all(IGeom,:)
143
           WRITE(*,*) 'Geom='
144
           WRITE(*,'(12(1X,F6.3))') Geom
145
           WRITE(*,*) 'DBG Egrad L165 GeomCart_old='
146
           DO I=1,Nat
147
              WRITE(*,'(3(1X,F12.8))') XyzGeomF(IGeom,:,I)
148
           END DO
149

    
150
           Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),Geom,XyzGeomF(IGeom,1,:), &
151
          XyzGeomF(IGeom,2,:),XyzGeomF(IGeom,3,:),x,y,z,XPrim,XPrimRef)
152
           if (debug) Print *, 'I am after Call ConvertBakerInternal_cart; L111, Egrad.f90'
153
           DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart()
154
              BTransInvF(IGeom,I,:) = BTransInv_local(I,:)
155
           END DO
156
           ! GeomCart is actually XyzCoord and has INTENT(OUT) in this subroutine.
157
           GeomCart(:,1) = x(:)
158
           GeomCart(:,2) = y(:)
159
           GeomCart(:,3) = z(:)
160
        END IF ! matches IF (IOpt .EQ. 0) THEN
161
     END IF
162
  CASE ('MIXED')
163
!  write(*,*) "Coucou 4-mixed"
164
     Call Mixed2Cart(Nat,indzmat,geom,GeomCart)
165
  CASE ('CART')
166

    
167
     GeomCart=reshape(Geom,(/Nat,3/))
168
     if (debug) THEN
169
        WRITE(*,*) "Coord=CART... in egrad"
170
        DO Iat=1,Nat
171
           WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),Geom(3*Iat-2:3*Iat)
172
        END DO
173
     END IF
174
  CASE ('HYBRID')
175
!  write(*,*) "Coucou 4-hybrid"
176
        GeomCart=reshape(Geom,(/Nat,3/))
177
  CASE DEFAULT
178
     WRITE (*,*) "Coord=",TRIM(COORD)," not recognized. Stop"
179
     STOP
180
  END SELECT
181
  !WRITE(*,*) Nat
182
  !WRITE(*,*) 'GeomCart in Egrad_baker'
183
  !DO I=1,Nat
184
  !  WRITE(*,'(1X,A,3(1X,F12.8))') AtName(I),GeomCart(I,1:3)
185
  !END DO
186

    
187
  SELECT CASE (Prog)
188
  CASE ('GAUSSIAN')
189
! we call the Gaussian routine.
190
! it will return the energy and the gradient in cartesian coordinates.
191
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT)
192
     Call egrad_gaussian(E,GeomCart,GradCart)
193
  CASE ('MOPAC')
194
! we call the Mopac routine.
195
! it will return the energy and the gradient in cartesian coordinates.
196
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT)
197
     Call egrad_mopac(E,GeomCart,GradCart)
198
  CASE ('VASP')
199
     Call egrad_vasp(E,Geomcart,GradCart)
200
  CASE ('SIESTA')
201
     Call egrad_siesta(E,Geomcart,GradCart)
202
  CASE ('TURBOMOLE')
203
     Call egrad_turbomole(E,Geomcart,GradCart)
204
  CASE ('EXT')
205
     Call egrad_ext(E,Geomcart,GradCart)
206
  CASE ('TEST')
207
     Call egrad_test(Nat,E,Geomcart,GradCart)
208
  CASE ('TEST2D')
209
     Call egrad_test_2D(Nat,E,Geomcart,GradCart)
210
  CASE ('CHAMFRE')
211
     Call egrad_chamfre(Nat,E,Geomcart,GradCart)
212
  CASE ('LEPS')
213
     Call egrad_LEPS(Nat,E,Geomcart,GradCart)
214
  END SELECT
215
  if (debug) THEN
216
     WRITE(*,*) 'DBG EGRAD, GradCart read'
217
     DO I=1,Nat
218
        Iat=I
219
        if (renum) Iat=orderInv(I)
220
        WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GradCart(3*I-2:3*I)
221
     END DO
222
  END IF
223

    
224
!  If (PROG=="VASP") STOP
225

    
226

    
227
  ! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid
228
  !  but we have to convert it into Zmat if COORD=Zmat
229
  SELECT CASE (COORD)
230
     CASE ("ZMAT")
231
        ALLOCATE(GradTmp(3*Nat))
232
        if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int"
233
        CALL CalcBMat_int (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome)
234

    
235
        if (debug) WRITE(*,*) "DBG EGRAD Calling ConvGrad_Cart2Int"
236
        CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
237

    
238
        if (debug) WRITE(*,*) "DBG EGRAD Storing Grad"
239
     Grad(1)=GradTmp(4)
240
     Grad(2)=GradTmp(7)
241
! We might have troubles whith angles that are not in the [0:pi] range because,
242
! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian...
243
! so that the derivative is not good, and a multiplicative factor should be added...
244
!
245
! This in fact should be taken care of in B mat calculation... 
246
!
247
     IF ((Geom(3).LE.0).OR.(Geom(3).GE.Pi)) THEN
248
        Grad(3)=-1.0d0*GradTmp(8)
249
     ELSE
250
        Grad(3)=GradTmp(8)
251
     END IF
252
     Idx=4
253
     DO I=4,Nat
254
        Grad(Idx)=GradTmp(3*i-2)
255
        Grad(Idx+1)=GradTmp(3*i-1)
256
        Grad(Idx+2)=GradTmp(3*i)
257
        Idx=Idx+3
258
     END DO
259
     DEALLOCATE(GradTmp)
260
   CASE ('BAKER')
261
     ! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid
262
     ! but we have to convert it into internal coordinates if COORD=Baker	 
263
!!!!!!!!!!!!!!!!!!!!
264
     !
265
     !   PFL Debug
266
     !
267
!!!!!!!!!!!!!!!!!!!!
268
     if (debugPFL) THEN
269
        WRITE(*,*) "DBG PFL Egrad_Baker... computing dzdc"
270
     if (.not.ALLOCATED(indzmat)) THEN
271
        ALLOCATE(indzmat(Nat,5))
272
     IndZmat=0
273
     Indzmat(1,1)=1
274
     IndZmat(2,1)=2
275
     IndZmat(2,2)=1
276
     IndZmat(3,1)=3
277
     IndZmat(3,2)=2
278
     IndZmat(3,3)=1
279
     IndZmat(4,1)=4
280
     IndZmat(4,2)=3
281
     IndZmat(4,3)=2
282
     IndZmat(4,4)=1
283
     END IF
284
     IF (.NOT.ALLOCATED(DzDc)) THEN
285
        ALLOCATE(DzDc(3,nat,3,Nat))
286
     END IF
287
     if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat))                                            
288
     DO I=1,Nat
289
     GeomTmp(1,I)=GeomCart(OrderInv(I),1)
290
     GeomTmp(2,I)=GeomCart(OrderInv(i),2)
291
     GeomTmp(3,I)=GeomCart(OrderInv(i),3)
292
     END DO
293

    
294
     DzDc=0.d0
295
     CALL CalcBMat_int (nat, GeomTmp, indzmat, dzdc, massat,atome)     
296

    
297
     END IF ! if (debugPFL)
298
!!!!!!!!!!!!!!!!!!!!!!!!!
299
     ! Debugging PFL
300
!!!!!!!!!!!!!!!!!!!!!!!!
301

    
302
     WRITE(*,*) "BTransInv_local Trans that is used originally"
303
     DO J=1,3*Nat
304
        WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J)
305
     END DO
306

    
307
     WRITE(*,*) "Calculating actual values using GeomCart"
308
     if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat))
309
     GeomTmp(1,:)=GeomCart(1:Nat,1)
310
     GeomTmp(2,:)=GeomCart(1:Nat,2)
311
     GeomTmp(3,:)=GeomCart(1:Nat,3)
312

    
313
     ! Computing B^prim:
314
     FAloBBT=.NOT.ALLOCATED(BBT)
315
     FAloBBTInv=.NOT.ALLOCATED(BBT_inv)
316
     FAloBPrimT=.NOT.ALLOCATED(BprimT)
317
     if (FAloBPrimT) ALLOCATE(Bprimt(3*Nat,NPrim))
318
     if (FAloBBT) ALLOCATE(BBt(NCoord,NCoord))
319
     if (FAloBBTInv) ALLOCATE(BBt_inv(NCoord,NCoord))
320
     BprimT=0.d0
321
     ScanCoord=>Coordinate
322
     I=0
323
     DO WHILE (Associated(ScanCoord%next))
324
        I=I+1
325
        SELECT CASE (ScanCoord%Type)
326
        CASE ('BOND') 
327
           CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,GeomTmp,BprimT(1,I))
328
        CASE ('ANGLE')
329
           CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,GeomTmp,BprimT(1,I))
330
        CASE ('DIHEDRAL')
331
           CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,GeomTmp,BprimT(1,I))
332
        END SELECT
333
        ScanCoord => ScanCoord%next
334
     END DO
335

    
336
     if (debug) THEN
337
        WRITE(*,*) "Bprim "
338
        DO J=1,3*Nat
339
           WRITE(*,'(50(1X,F12.6))') BprimT(J,:)
340
        END DO
341
     END IF
342

    
343
     BMat_BakerT = 0.d0
344
     DO I=1,NCoord
345
        DO J=1,NPrim
346
           !BprimT is transpose of B^prim.
347
           !B = UMat^T B^prim, B^T = (B^prim)^T UMat
348
           BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I)
349
        END DO
350
     END DO
351

    
352
     BBT = 0.d0
353
     DO I=1,NCoord
354
        DO J=1,3*Nat   ! BBT(:,I) forms BB^T
355
           BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I)
356
        END DO
357
     END DO
358

    
359
     BBT_inv = 0.d0
360

    
361
     Call GenInv(NCoord,BBT,BBT_inv,NCoord)
362

    
363
!     ALLOCATE(V(NCoord,NCoord))
364
!     tol = 1e-10
365
!     ! BBT is destroyed by GINVSE
366
!     Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT)
367
!     DEALLOCATE(V)
368
!     IF(Fail) Then
369
!        WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90"
370
!        STOP
371
!     END IF
372

    
373
     !Print *, 'BBT_inv='
374
     DO J=1,NCoord
375
        ! WRITE(*,'(50(1X,F10.2))') BBT_inv(:,J)
376
        ! Print *, BBT_inv(:,J)
377
     END DO
378
     !Print *, 'End of BBT_inv'
379

    
380
     ! Calculation of (B^T)^-1 = (BB^T)^-1B:
381
     BTransInv_local = 0.d0
382
     DO I=1,3*Nat
383
        DO J=1,NCoord
384
           BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J)
385
        END DO
386
     END DO
387

    
388
     Print *, 'BMat_BakerT='
389
     DO J=1,3*Nat
390
        WRITE(*,'(50(1X,F12.8))') BMat_BakerT(:,J)
391
     END DO
392

    
393
     WRITE(*,*) "BTransInv_local Trans corresponding to GeomCart"
394
     DO J=1,3*Nat
395
        WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J)
396
     END DO
397

    
398
     if (FAloBPrimT) DEALLOCATE(Bprimt)
399
     if (FAloBBT) DEALLOCATE(BBt)
400
     if (FAloBBTInv) DEALLOCATE(BBt_inv)
401

    
402
     Grad=0.d0
403
     DO I=1, 3*Nat
404
        ! Here Grad is gradient in Baker coordinates. GradCart(3*Nat)
405
        Grad(:) = Grad(:) + BTransInv_local(:,I)*GradCart(I) !BTransInv_local = BTransInv in Opt_Geom.f90
406
     END DO
407
     !WRITE(*,*) 'In Egrad.f90, L263, Gradient:' 
408
     !WRITE(*,'(12(1X,F8.4))') Grad(1:NCoord)
409

    
410

    
411
     CASE ('MIXED')
412
        ALLOCATE(GradTmp(3*Nat))
413
        if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed"
414
        CALL CalcBMat_mixed (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome)
415

    
416
        if (Debug) THEN
417
           WRITE(*,*) "DzDc"
418
           DO I=1,Nat
419
              DO J=1,3
420
                 WRITE(*,'(50(1X,F12.6))') dzdc(J,I,:,:)
421
              END DO
422
           END DO
423
        END IF
424
       
425

    
426
        CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
427
        DO I=1,Nat
428
           WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I-2:3*I)
429
        END DO
430
        SELECT CASE (NCART)
431
           CASE (1)
432
              Grad(1:3)=GradTmp(1:3)
433
              Grad(4)=GradTmp(4)
434
              Grad(5)=GradTmp(7)
435
              Grad(6)=GradTmp(8)
436
              Idx=7
437
              IBeg=4
438
           CASE(2)
439
              Grad(1:3)=GradTmp(1:3)
440
              Grad(4:6)=GradTmp(4:6)
441
              Grad(7)=GradTmp(7)
442
              Grad(8)=GradTmp(8)
443
              Idx=9
444
               IBeg=4
445
            CASE DEFAULT
446
               Idx=1
447
               IBeg=1
448
            END SELECT
449
     DO I=IBeg,Nat
450
        Grad(Idx)=GradTmp(3*i-2)
451
        Grad(Idx+1)=GradTmp(3*i-1)
452
        Grad(Idx+2)=GradTmp(3*i)
453
        Idx=Idx+3
454
     END DO
455
     DEALLOCATE(GradTmp)
456
  CASE ("CART","HYBRID")
457
! PFL 2013 Feb 
458
! We could think of just doing:
459
!     Grad=GradCart
460
! however, this does not work as implicitly we have :
461
! Grad(Nat,3)  BUT GradCar(3,Nat)
462
! So we have to use reshape :-(
463
     Grad=reshape(reshape(GradCart,(/3,Nat/),ORDER=(/2,1/)),(/3*Nat/))
464
  CASE DEFAULT
465
     WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in Egradients. STOP"
466
     STOP
467
  END SELECT
468

    
469
  DEALLOCATE(GradCart)
470
  DEALLOCATE(x,y,z)
471

    
472
  IF (Debug) WRITE(*,*) 'DBG EGRAD grad',grad(1:NCoord)
473
  if (debug) Call Header("Egrad Over")
474

    
475
  !999 CONTINUE
476
  !if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!'
477
  !STOP
478
  ! ======================================================================
479
end subroutine Egrad
480