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