Statistiques
| Révision :

root / src / ConvertBakerInternal_cart.f90 @ 8

Historique | Voir | Annoter | Télécharger (23,95 ko)

1
SUBROUTINE ConvertBakerInternal_cart(IntCoord_k,IntCoord,x_k,y_k,z_k,x,y,z,XPrim,XPrimRef)
2

    
3
!================================================================
4
! Converts the positions in Baker Coordinates to cartesian coordinates
5
!================================================================
6
!
7
! Input:
8
!    Incoord_k(NCoord) REAL: Starting internal coordinates
9
!    IntCoord(NCoord) REAL: Coordinate to convert to cartesian coordinates
10
!    x_k,y_k,z_k(Nat) REAL: Starting cartesian geometry
11
!    XPrimRef(NPrim) REAL: Reference values for Primitives coordinates (mainly for dihedral)
12
!
13
! Output:
14
!    x,y,z(Nat) REAL: Cartesian coordinates corresponding to IntCoord
15
!    XPrim(NPrim) REAL: Primitives for the cartesian coordinates. (This is a by-product)
16
!
17
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
18

    
19
! BTransInv,UMat should be passed as arguments of this subroutine.
20

    
21

    
22
  ! IdxGeom: geometry index.
23
  ! Internal coordinates IntCoord(NCoord) is to be converted
24
  ! into Cartesian coordinates x(Nat),y(Nat),z(Nat).
25

    
26
  use Path_module, only : Nat,AtName,BMat_BakerT,NPrim,&
27
       BBT,BBT_inv,BprimT,ScanCoord,Coordinate,&
28
       BTransInv_local,Xprimitive_t,&
29
       NCoord,UMat_local
30
  ! UMat_local is allocated in Calc_baker.f90
31
  ! XyzGeomI(NGeomI,3,Nat),BMat_BakerT(3*Nat,NCoord),a0=0.529177249d0
32

    
33
  Use Io_Module
34
  use VarTypes
35
  IMPLICIT NONE
36

    
37
  REAL(KREAL), INTENT(OUT) :: x(Nat),y(Nat),z(Nat)
38
  REAL(KREAL), INTENT(IN) :: x_k(Nat),y_k(Nat),z_k(Nat)
39
  REAL(KREAL), INTENT(IN) :: IntCoord(NCoord)
40
  REAL(KREAL), INTENT(INOUT) ::  IntCoord_k(NCoord)
41
  REAL(KREAL), INTENT(IN) :: XPrimRef(NPrim)
42
  REAL(KREAL), INTENT(OUT) :: XPrim(NPrim)
43

    
44

    
45
  ! IdxGeom: geometry index.
46
  Integer(KINT) :: I, J, Counter
47
  REAL(KREAL) :: normDeltaX_int
48

    
49
  REAL(KREAL), ALLOCATABLE :: DeltaX_int(:)  !DeltaX_int(NPrim,3*Nat-6), 3*Nat-6??
50
  REAL(KREAL), ALLOCATABLE :: DeltaX_cart(:) !DeltaX_cart(3*Nat)
51
  REAL(KREAL), ALLOCATABLE :: Geom(:,:) !(3,Nat)
52
  REAL(KREAL), ALLOCATABLE :: X_cart_k(:,:) ! (3,Nat)
53
  REAL(KREAL), ALLOCATABLE :: GMat(:,:) !(NPrim,NPrim)
54
  REAL(KREAL), ALLOCATABLE :: EigVec(:,:), EigVal(:) ! EigVec(NPrim,NPrim)
55
  REAL(KREAL), ALLOCATABLE :: UMat_tmp(:,:)
56

    
57
  REAL(KREAL) :: Angle, angle_d
58
  REAL(KREAL), PARAMETER :: Crit=1e-08
59

    
60
  EXTERNAL Angle, angle_d
61
  LOGICAL :: debug, FAIL
62

    
63
  INTEGER(KINT), PARAMETER :: NbItMax=50
64

    
65

    
66
  INTERFACE
67
     function valid(string) result (isValid)
68
       CHARACTER(*), intent(in) :: string
69
       logical                  :: isValid
70
     END function VALID
71

    
72

    
73

    
74
     SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef)
75
       !     
76
       ! This subroutine uses the description of a list of Coordinates
77
       ! to compute the values of the coordinates for a given geometry.
78
       !
79
!!!!!!!!!!
80
       ! Input:
81
       !   Na: INTEGER, Number of atoms
82
       !  x,y,z(Na): REAL, cartesian coordinates of the considered geometry
83
       ! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes
84
       ! NPrim, INTEGER: Number of coordinates to compute
85
       !
86
       ! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for
87
       ! a former geometry. Useful for Dihedral angles evolution...
88

    
89
!!!!!!!!!!!
90
       ! Output:
91
       ! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates
92
       !
93
!!!!!!!!!
94

    
95
       Use VarTypes
96
       Use Io_module
97
       Use Path_module, only : pi
98

    
99
       IMPLICIT NONE
100

    
101
       Type (ListCoord), POINTER :: Coordinate
102
       INTEGER(KINT), INTENT(IN) :: Nat,NPrim
103
       REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat)
104
       REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) 
105
       REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim)
106

    
107
     END SUBROUTINE CALC_XPRIM
108
  END INTERFACE
109

    
110

    
111
  debug=valid('ConvertBakerInternal_cart')
112
  if (debug) WRITE(*,*) '=======Entering ConvertBakerInternal_cart======='
113

    
114
  FAIL = .FALSE.     
115

    
116
!!!!!!
117
  ! Allocation phase
118

    
119
  ALLOCATE(DeltaX_int(NCoord),DeltaX_cart(3*Nat))
120
  ALLOCATE(UMat_tmp(NPrim,NCoord))
121
  ALLOCATE(X_cart_k(3,Nat))
122
  ALLOCATE(Geom(3,Nat),BprimT(3*Nat,NPrim))
123
  ALLOCATE(BBT(NCoord,NCoord))
124
  ALLOCATE(BBT_inv(NCoord,NCoord))
125
  ALLOCATE(Gmat(NPrim,NPrim))
126
  ALLOCATE(EigVal(NPrim),EigVec(NPrim,NPrim))
127

    
128

    
129
  if (debug) THEN
130
     WRITE(*,*) "Checking purposes...."
131
     WRITE(*,*) "Trying ot convert Xint initial (IntCoord)"
132
     WRITE(*,'(50(1X,F12.8))') IntCoord(:)
133
     WRITE(*,*) "BTransInv_local"
134
     DO J=1,3*Nat
135
        WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J)
136
     END DO
137
     WRITE(*,*) "Xint initial (IntCoord_k)"
138
     WRITE(*,'(50(1X,F12.8))') IntCoord_k(:)
139
     WRITE(*,*) "Cartesian coordinates"
140
     WRITE(*,*) Nat
141
     WRITE(*,*) 'GeomCart in ConvertBakerInternal_cart'
142
     DO I=1,Nat
143
        WRITE(*,'(1X,A,I4,3(1X,F12.4))') AtName(I),I,X_k(I),Y_k(i),Z_k(i)
144
     END DO
145
     WRITE(*,*) "Calculating actual values using x_k,y_k,z_k"
146

    
147
!     WRITE(*,*) "First, with Calc_XPrim"
148
!     WRITE(*,*) "Coordinate:",associated(Coordinate)
149
     Call Calc_Xprim(nat,x_k,y_k,z_k,Coordinate,NPrim,XPrimitive_t,XPrimRef)
150

    
151
!      I=0 ! index for the NPrim (NPrim is the number of ! primitive coordinates).
152
!      ScanCoord=>Coordinate
153
!      DO WHILE (Associated(ScanCoord%next))
154
!         I=I+1
155
!         SELECT CASE (ScanCoord%Type)
156
!         CASE ('BOND')
157
!            Print *, I, ':',  ScanCoord%At1, '-', ScanCoord%At2, '=', Xprimitive_t(I)
158
!         CASE ('ANGLE')    
159
!            Print *, I, ':', ScanCoord%At1, '-', ScanCoord%At2,'-',ScanCoord%At3, '=', Xprimitive_t(I)*180./Pi
160
!         CASE ('DIHEDRAL')
161
!            Print *, I, ':',  ScanCoord%At1, '-',ScanCoord%At2,'-',ScanCoord%At3,'-',&
162
!            ScanCoord%At4,'=',Xprimitive_t(I)*180./Pi
163
!         END SELECT
164
!         ScanCoord => ScanCoord%next
165
!      END DO ! matches DO WHILE (Associated(ScanCoord%next))
166

    
167
!      WRITE(*,*) "Second with normal proc"
168
!      I=0 ! index for the NPrim (NPrim is the number of ! primitive coordinates).
169
!      Xprimitive_t=0.d0
170
!      ScanCoord=>Coordinate
171
!      DO WHILE (Associated(ScanCoord%next))
172
!         I=I+1
173
!         SELECT CASE (ScanCoord%Type)
174
!         CASE ('BOND')
175
!            Call vecteur(ScanCoord%At2,ScanCoord%At1,x_k,y_k,z_k,vx2,vy2,vz2,Norm2)
176
!            Xprimitive_t(I) = Norm2
177
!            Print *, I, ':',  ScanCoord%At1, '-', ScanCoord%At2, '=', Xprimitive_t(I)
178
!         CASE ('ANGLE')    
179
!            Call vecteur(ScanCoord%At2,ScanCoord%At3,x_k,y_k,z_k,vx1,vy1,vz1,Norm1)
180
!            Call vecteur(ScanCoord%At2,ScanCoord%At1,x_k,y_k,z_k,vx2,vy2,vz2,Norm2)
181
!            Xprimitive_t(I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
182
!            Print *, I, ':', ScanCoord%At1, '-', ScanCoord%At2,'-',ScanCoord%At3, '=', Xprimitive_t(I)*180./Pi
183
!         CASE ('DIHEDRAL')
184
!            Call vecteur(ScanCoord%At2,ScanCoord%At1,x_k,y_k,z_k,vx1,vy1,vz1,Norm1)
185
!            Call vecteur(ScanCoord%At2,ScanCoord%At3,x_k,y_k,z_k,vx2,vy2,vz2,Norm2)
186
!            Call vecteur(ScanCoord%At3,ScanCoord%At4,x_k,y_k,z_k,vx3,vy3,vz3,Norm3)
187
!            Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5)
188
!            Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
189
! !!! PFL 28 Feb 2008 Test to be sure that angle does not change by more than Pi
190
!            !           Xprimitive_t(I)=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5,vx2,vy2,vz2,norm2)*Pi/180.
191
!            DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
192
!                 vx2,vy2,vz2,norm2)
193
!            Xprimitive_t(I) = DiheTmp*Pi/180.
194
!            ! We treat large dihedral angles differently as +180=-180 mathematically and physically
195
!            ! but this causes lots of troubles when converting baker to cart.
196
!            ! So we ensure that large dihedral angles always have the same sign
197
!            if (abs(ScanCoord%SignDihedral).NE.1) THEN
198
!               ScanCoord%SignDihedral=Int(Sign(1.,DiheTmp))
199
!            ELSE
200
!               If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN
201
!                  Xprimitive_t(I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi
202
!               END IF
203
!            END IF
204

    
205
!            IF (abs(Xprimitive_t(I)-XPrimRef(I)).GE.Pi) THEN
206
!               if ((Xprimitive_t(I)-XPrimRef(I)).GE.Pi) THEN
207
!                  Xprimitive_t(I)=Xprimitive_t(I)-2.d0*Pi
208
!               ELSE
209
!                  Xprimitive_t(I)=Xprimitive_t(I)+2.d0*Pi
210
!               END IF
211
!            END IF
212

    
213
!            Print *, I, ':',  ScanCoord%At1, '-',ScanCoord%At2,'-',ScanCoord%At3,'-',&
214
!            ScanCoord%At4,'=',Xprimitive_t(I)*180./Pi
215
!         END SELECT
216
!         ScanCoord => ScanCoord%next
217
!      END DO ! matches DO WHILE (Associated(ScanCoord%next))
218

    
219
     IntCoord_k=0.d0
220
     DO I=1, NPrim
221
        ! Transpose of UMat is needed below, that is why UMat(I,:).
222
        IntCoord_k(:)=IntCoord_k(:)+UMat_local(I,:)*Xprimitive_t(I)     
223
     END DO
224
     WRITE(*,*) "Xint (intCoord_k) cooresponding to x_k,y_k,z_k"
225
     WRITE(*,'(50(1X,F12.8))') IntCoord_k(:)
226
     WRITE(*,*) "UMat_local"
227
     DO I=1, NPrim
228
        WRITE(*,'(50(1X,F12.8))') UMat_local(I,:)
229
     END DO
230

    
231
     !     Geom(1,:)=X_k(1:Nat)/a0
232
     !     Geom(2,:)=y_k(1:Nat)/a0
233
     !     Geom(3,:)=z_k(1:Nat)/a0
234

    
235
     Geom(1,:)=X_k(1:Nat)
236
     Geom(2,:)=y_k(1:Nat)
237
     Geom(3,:)=z_k(1:Nat)
238

    
239
     ! Computing B^prim:
240
     BprimT=0.d0
241
     ScanCoord=>Coordinate
242
     I=0
243
     DO WHILE (Associated(ScanCoord%next))
244
        I=I+1
245
        SELECT CASE (ScanCoord%Type)
246
        CASE ('BOND') 
247
           CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,Geom,BprimT(1,I))
248
        CASE ('ANGLE')
249
           CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,Geom,BprimT(1,I))
250
        CASE ('DIHEDRAL')
251
           CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I))
252
        END SELECT
253
        ScanCoord => ScanCoord%next
254
     END DO
255

    
256
     !     if (debug) THEN
257
     !        WRITE(*,*) "PFL DBG, I,NPrim,BPrimT",I,NPrim
258
     !        DO I=1,NPrim
259
     !           WRITE(*,'(50(1X,F12.6))') BPrimT(:,I)
260
     !        END DO
261
     !     END IF
262

    
263
     BMat_BakerT = 0.d0
264
     DO I=1,NCoord
265
        DO J=1,NPrim
266
           !BprimT is transpose of B^prim.
267
           !B = UMat^T B^prim, B^T = (B^prim)^T UMat
268
           BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I)
269
        END DO
270
     END DO
271

    
272
     BBT = 0.d0
273
     DO I=1,NCoord
274
        DO J=1,3*Nat   ! BBT(:,I) forms BB^T
275
           BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I)
276
        END DO
277
     END DO
278

    
279
     BBT_inv = 0.d0
280

    
281
     !Print *, 'NCoord=', NCoord
282
     !Print *, 'BBT='
283
     DO J=1,NCoord
284
        ! WRITE(*,'(50(1X,F12.6))') BBT(:,J)
285
     END DO
286
     !Print *, 'End of BBT'
287
     ! GenInv is in Mat_util.f90
288
     Call GenInv(NCoord,BBT,BBT_inv,NCoord)
289
!     ALLOCATE(V(NCoord,NCoord))
290
!     tol = 1e-12
291
!     ! BBT is destroyed by GINVSE
292
!     Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT)
293
!     DEALLOCATE(V)
294
!     IF(Fail) Then
295
!        WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90"
296
!        STOP
297
!     END IF
298

    
299
     !Print *, 'BBT_inv='
300
     DO J=1,NCoord
301
        ! WRITE(*,'(50(1X,F10.2))') BBT_inv(:,J)
302
        ! Print *, BBT_inv(:,J)
303
     END DO
304
     !Print *, 'End of BBT_inv'
305

    
306
     ! Calculation of (B^T)^-1 = (BB^T)^-1B:
307
     BTransInv_local = 0.d0
308
     DO I=1,3*Nat
309
        DO J=1,NCoord
310
           BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J)
311
        END DO
312
     END DO
313

    
314
     !Print *, 'BMat_BakerT='
315
     DO J=1,3*Nat
316
        !WRITE(*,'(50(1X,F12.8))') BMat_BakerT(:,J)
317
     END DO
318
     !Print *, 'End of BMat_BakerT'
319

    
320
     WRITE(*,*) "BTransInv_local Trans corresponding to x_k,y_k,z_k"
321
     DO J=1,3*Nat
322
        WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J)
323
     END DO
324

    
325

    
326
  END IF
327

    
328
  DeltaX_int(:) = IntCoord(:)-IntCoord_k(:)
329

    
330
  X_cart_k(1,:) = x_k(:)
331
  X_cart_k(2,:) = y_k(:)
332
  X_cart_k(3,:) = z_k(:)
333

    
334
  normDeltaX_int=0.d0 ! This is to be implemented.
335
  DO I=1, NCoord
336
     normDeltaX_int=normDeltaX_int+DeltaX_int(I)*DeltaX_int(I)
337
  END DO
338
  normDeltaX_int = sqrt(normDeltaX_int)
339
  ! I just need initial value of normDeltaX_int to be .GT. 1d-11,
340
  ! that is why it is initialized to 1.d0
341
  !normDeltaX_int = 1.d0
342

    
343
  if (debug) WRITE(*,*) "Entering the loop"
344
  Counter = 0
345
  DO While ((normDeltaX_int .GT. Crit) .AND. (Counter .LT. NbItMax))
346
     !DO While (normDeltaX_int .GT. 1d-6)
347
     Counter = Counter + 1       
348
     DeltaX_cart = 0.d0
349
!     if (normDeltaX_int.LE.1e-4) THEN
350
!        DeltaX_int=DeltaX_int*1e4
351
!     END IF 
352
     DO I=1,NCoord
353
        ! below BTransInv_local(I,:) is used because actually we need Transpose of BTransInv_local.
354
        DeltaX_cart(:)=DeltaX_cart(:)+BTransInv_local(I,:)*DeltaX_int(I)
355
     END DO
356
!     if (normDeltaX_int.LE.1e-4) THEN
357
!        DeltaX_int=DeltaX_int*1e-4
358
!     DeltaX_cart=DeltaX_cart*1e-4
359
!     END IF
360

    
361

    
362
           if (debug) THEN
363
              WRITE(*,*) "Info for iteration:",counter
364
              Print *, 'DeltaX_int='
365
              WRITE(*,'(50(1X,F12.8))') DeltaX_int
366
              Print *, 'DeltaX_cart,norm',sqrt(dot_product(DeltaX_cart,DeltaX_cart))
367
              DO J=1,Nat
368
              WRITE(*,'(1X,A4,3(1X,F15.11),3(1X,D25.15))') AtName(J),DeltaX_cart(3*J-2:3*J),DeltaX_cart(3*J-2:3*J)
369
              END DO
370
              Print *, 'BTransInv_local Trans='
371
              DO J=1,3*Nat
372
                 WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J)
373
              END DO
374
              Print *, 'DeltaX_int='
375
              WRITE(*,'(50(1X,F12.8))') DeltaX_int
376
           END IF
377

    
378
     ! Finite precision error correction step:
379
!     DO I=1,3*Nat
380
!        IF (DeltaX_cart(I) .NE. 0.d0) Then
381
!           IF(abs(DeltaX_cart(I)) .LT. 1d-10) Then
382
!              !Print *, 'Correcting DeltaX_cart(',I,')=', DeltaX_cart(I)
383
!              DeltaX_cart(I) = 0.d0
384
!           END IF
385
!        END IF
386
!     END DO
387

    
388
     ! new cartesian coordinates:
389
     DO I=1, Nat
390
        X_cart_k(1,I) = X_cart_k(1,I) + DeltaX_cart(3*I-2)
391
        X_cart_k(2,I) = X_cart_k(2,I) + DeltaX_cart(3*I-1)
392
        X_cart_k(3,I) = X_cart_k(3,I) + DeltaX_cart(3*I)
393
     END DO
394

    
395
     !     Geom(1,:)=X_cart_k(1,1:Nat)/a0
396
     !     Geom(2,:)=X_cart_k(2,1:Nat)/a0
397
     !     Geom(3,:)=X_cart_k(3,1:Nat)/a0
398
     Geom(1,:)=X_cart_k(1,1:Nat)
399
     Geom(2,:)=X_cart_k(2,1:Nat)
400
     Geom(3,:)=X_cart_k(3,1:Nat)
401

    
402
     if (debug) THEN
403
     WRITE(*,*) 'GeomCart used to compute BPrim'
404
     DO I=1,Nat
405
        WRITE(*,'(1X,A,I4,3(1X,F12.4))') AtName(I),I,Geom(1,I),Geom(2,i),Geom(3,i)
406
     END DO
407
     END IF
408

    
409
     ! Computing B^prim:
410
     BprimT=0.d0
411
     ScanCoord=>Coordinate
412
     I=0
413
     DO WHILE (Associated(ScanCoord%next))
414
        I=I+1
415
        SELECT CASE (ScanCoord%Type)
416
        CASE ('BOND') 
417
           CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,Geom,BprimT(1,I))
418
        CASE ('ANGLE')
419
           CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,Geom,BprimT(1,I))
420
        CASE ('DIHEDRAL')
421
           CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I))
422
        END SELECT
423
        ScanCoord => ScanCoord%next
424
     END DO
425

    
426
     if (debug) THEN
427
        WRITE(*,*) "Bprim "
428
        DO J=1,3*Nat
429
           WRITE(*,'(50(1X,F12.6))') BprimT(J,:)
430
        END DO
431
     END IF
432

    
433

    
434
     BMat_BakerT = 0.d0
435
     DO I=1,NCoord
436
        DO J=1,NPrim
437
           !BprimT is transpose of B^prim.
438
           !B = UMat^T B^prim, B^T = (B^prim)^T UMat
439
           BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I)
440
        END DO
441
     END DO
442

    
443
     ! ADDED FROM HERE:
444
     ! We now compute G=B(BT) matrix
445
     !IF (IOptt .EQ. 5000) Then
446
     !GMat=0.d0
447
     !DO I=1,NPrim
448
     !  DO J=1,3*Nat
449
     !    GMat(:,I)=Gmat(:,I)+BprimT(J,:)*BprimT(J,I) !*1.d0/mass(atome(int(K/3.d0)))
450
     !END DO
451
     !END DO
452

    
453
     ! We diagonalize G
454
     !Call Jacobi(GMat,NPrim,EigVal,EigVec,NPrim)
455
     !Call Trie(NPrim,EigVal,EigVec,NPrim)
456

    
457
     ! We construct the new DzDc(b=baker) matrix
458
     ! UMat is nonredundant vector set, i.e. set of eigenvectors of
459
     !  BB^T corresponding to eigenvalues > zero.
460

    
461
     !UMat=0.d0
462
     !UMat_tmp=0.d0
463
     !BMat_BakerT = 0.d0
464
     !J=0
465
     !DO I=1,NPrim
466
     !  IF (abs(eigval(I))>=1e-6) THEN
467
     !   J=J+1
468
     !  DO K=1,NPrim
469
     ! BprimT is transpose of B^prim.
470
     ! B = UMat^T B^prim, B^T = (B^prim)^T UMat
471
     !BMat_BakerT(:,J)=BMat_BakerT(:,J)+BprimT(:,K)*Eigvec(K,I)
472
     ! END DO
473
     !IF(J .GT. 3*Nat-6) THEN
474
     !WRITE(*,*) 'Number of vectors in Eigvec with eigval .GT. &
475
     !    1e-6 (=UMat) (=' ,J,') exceeded 3*Nat-6=',3*Nat-6, &
476
     !    'Stopping the calculation.'
477
     !   STOP
478
     !  END IF
479
     ! UMat_tmp(:,J) =  Eigvec(:,I)
480
     ! END IF
481
     !END DO
482

    
483
     ! THIS IS TO BE CORRECTED, A special case for debugging C2H3F case:
484
     !J=0
485
     !DO I=1, 3*Nat-6
486
     ! IF ((I .NE. 6) .AND. (I .NE. 7) .AND. (I .NE. 9)) Then
487
     !  J=J+1
488
     ! Print *, 'J=', J, ' I=', I
489
     !UMat(:,J) = UMat_tmp(:,I)
490
     !     END IF
491
     ! END DO
492
     !BMat_BakerT=0.d0
493
     !    DO I=1,NCoord
494
     !      DO J=1,NPrim
495
     ! B = UMat^T B^prim, B^T = (B^prim)^T UMat
496
     !       BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat(J,I)
497
     !    END DO
498
     !END DO
499
     ! TO BE CORRECTED, ENDS HERE.
500

    
501
     !END IF
502
     ! END of the new changes.
503

    
504
     BBT = 0.d0
505
     DO I=1,NCoord
506
        DO J=1,3*Nat   ! BBT(:,I) forms BB^T
507
           BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I)
508
        END DO
509
     END DO
510

    
511
     BBT_inv = 0.d0
512

    
513
     !Print *, 'NCoord=', NCoord
514
     !Print *, 'BBT='
515
     DO J=1,NCoord
516
        ! WRITE(*,'(50(1X,F12.6))') BBT(:,J)
517
     END DO
518
     !Print *, 'End of BBT'
519
     ! GenInv is in Mat_util.f90
520
     Call GenInv(NCoord,BBT,BBT_inv,NCoord)
521
     !     if (debug) THEN
522
     !        WRITE(*,*) 'BBT_inv by GenInv'
523
     !        DO J=1,NCoord
524
     !           WRITE(*,'(50(1X,F12.6))') BBT_inv(:,J)
525
     !        END DO
526
     !     END IF
527

    
528
!    ALLOCATE(V(NCoord,NCoord))
529
!    tol = 1e-12
530
!    ! BBT is destroyed by GINVSE
531
!    Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT)
532
!    DEALLOCATE(V)
533
!    IF(Fail) Then
534
!       WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90"
535
!       STOP
536
!    END IF
537

    
538
     !     if (debug) THEN
539
     !        WRITE(*,*) 'BBT_inv by Genvse'
540
     !        DO J=1,NCoord
541
     !           WRITE(*,'(50(1X,F12.6))') BBT_inv(:,J)
542
     !        END DO
543
     !     END IF
544

    
545
     !Print *, 'End of BBT_inv'
546

    
547
     ! Calculation of (B^T)^-1 = (BB^T)^-1B:
548
     BTransInv_local = 0.d0
549
     DO I=1,3*Nat
550
        DO J=1,NCoord
551
           BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J)
552
        END DO
553
     END DO
554

    
555
     !     if (debug) THEN
556
     !        Print *, 'BMat_BakerT='
557
     !        DO J=1,3*Nat
558
     !           WRITE(*,'(50(1X,F12.6))') BMat_BakerT(:,J)
559
     !        END DO
560
     !        Print *, 'End of BMat_BakerT'
561
     !     END IF
562

    
563
     !      if (debug) THEN
564
     !         Print *, 'BTransInv_local Trans='
565
     !         DO J=1,3*Nat
566
     !            WRITE(*,'(50(1X,F16.6))') BTransInv_local(:,J)
567
     !         END DO
568
     !         Print *, 'End of BTransInv_local Trans'
569
     !      END IF
570

    
571
     ! New cartesian cooordinates:
572
     x(:) = X_cart_k(1,:)
573
     y(:) = X_cart_k(2,:)
574
     z(:) = X_cart_k(3,:)
575

    
576

    
577
     Call Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive_t,XPrimRef)
578

    
579
!      I=0 ! index for the NPrim (NPrim is the number of 
580
!      ! primitive coordinates).
581
!      Xprimitive_t=0.d0
582
!      ScanCoord=>Coordinate
583
!      DO WHILE (Associated(ScanCoord%next))
584
!         I=I+1
585
!         SELECT CASE (ScanCoord%Type)
586
!         CASE ('BOND')
587
!            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)
588
!            Xprimitive_t(I) = Norm2
589
!            !Print *, I, ':',  ScanCoord%At1, '-', ScanCoord%At2, '=', Xprimitive_t(I)
590
!         CASE ('ANGLE')    
591
!            Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx1,vy1,vz1,Norm1)
592
!            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)
593
!            Xprimitive_t(I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
594
!            !Print *, I, ':', ScanCoord%At1, '-', ScanCoord%At2,'-',ScanCoord%At3, '=', Xprimitive_t(I)*180./Pi
595
!         CASE ('DIHEDRAL')
596
!            Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2)
597
!            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1)
598
!            Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3)
599
!            Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5)
600
!            Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
601
! !!! PFL 28 Feb 2008 Test to be sure that angle does not change by more than Pi
602
!            !           Xprimitive_t(I)=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5,vx2,vy2,vz2,norm2)*Pi/180.
603
!            DiheTmp=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5,vx2,vy2,vz2,norm2)
604
!            Xprimitive_t(I) = DiheTmp*Pi/180.
605
!            ! We treat large dihedral angles differently as +180=-180 mathematically and physically
606
!            ! but this causes lots of troubles when converting baker to cart.
607
!            ! So we ensure that large dihedral angles always have the same sign
608
!            if (abs(ScanCoord%SignDihedral).NE.1) THEN
609
!               ScanCoord%SignDihedral=Int(Sign(1.d0,DiheTmp))
610
!            ELSE
611
!               If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN
612
!                  Xprimitive_t(I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi
613
!               END IF
614
!            END IF
615

    
616

    
617

    
618
!            !Print *, I, ':',  ScanCoord%At1, '-',ScanCoord%At2,'-',ScanCoord%At3,'-',&
619
!            !ScanCoord%At4,'=',Xprimitive_t(I)*180./Pi
620
!         END SELECT
621
!         ScanCoord => ScanCoord%next
622
!      END DO ! matches DO WHILE (Associated(ScanCoord%next))
623

    
624
     !     if (debug) THEN
625
     !        WRITE(*,*) "Xprimitive_t"
626
     !        WRITE(*,'(50(1X,F12.6))') Xprimitive_t
627
     !     END IF
628

    
629
     IntCoord_k=0.d0
630
     DO I=1, NPrim
631
        ! Transpose of UMat is needed below, that is why UMat(I,:).
632
        IntCoord_k(:)=IntCoord_k(:)+UMat_local(I,:)*Xprimitive_t(I)     
633
     END DO
634

    
635
     if (debug) THEN
636
        WRITE(*,*) "New Xint (IntCoord_k)"
637
        WRITE(*,'(50(1X,F12.8))') IntCoord_k
638
        WRITE(*,*) "Target (IntCoord)"
639
        WRITE(*,'(50(1X,F12.8))') IntCoord
640

    
641
     END IF
642

    
643

    
644
     ! Computing DeltaX_int
645

    
646
     DeltaX_int(:) = IntCoord(:)-IntCoord_k(:)
647

    
648

    
649

    
650

    
651
     ! norm2 of DeltaX_int is being calculated here:
652
     normDeltaX_int=0.d0
653
     DO I=1, NCoord
654
        normDeltaX_int=normDeltaX_int+DeltaX_int(I)*DeltaX_int(I)
655
     END DO
656
     normDeltaX_int = sqrt(normDeltaX_int)
657

    
658
          if (debug) THEN
659
             WRITE(*,*) "New Delta_Xint (deltaX_int)"
660
             WRITE(*,'(50(1X,F12.6))') DeltaX_int
661
             WRITE(*,*) "Norm:",normDeltaX_int
662
          END IF
663

    
664
     !IF (Counter .GE. 400) THEN
665
     !Print *, 'Counter=', Counter, 'normDeltaX_int=', normDeltaX_int
666
     !Print *, 'DeltaX_int='
667
     !WRITE(*,'(50(1X,F8.2))') DeltaX_int
668
     !END IF
669

    
670
     IF (Mod(Counter,1).EQ.0) THEN
671
        !WRITE(*,*) Nat
672
        !WRITE(*,*) 'GeomCart in ConvertBakerInternal_cart'
673
        DO I=1,Nat
674
           !  WRITE(*,'(1X,A,I4,3(1X,F12.4))') AtName(I),I,X_Cart_k(1:3,I)
675
        END DO
676
     END IF
677

    
678
     !call two_norm(DeltaX_int,normDeltaX_int,3*Nat)
679
  END DO !matches DO While (normDeltaX_int .GT. 1d-11)
680

    
681
  IF (debug) THEN
682
  WRITE(*,*) 'GeomCart in ConvertBakerInternal_cart'
683
  DO I=1,Nat
684
     WRITE(*,'(1X,A,I4,3(1X,F12.4))') AtName(I),I,X_Cart_k(1:3,I)
685
  END DO
686
  END IF
687

    
688
  if (debug) THEN
689
     WRITE(*,*) "Bmat_bakerT"
690
     DO J=1,NCoord
691
        WRITE(*,'(50(1X,F12.6))') BMat_BakerT(:,J)
692
     END DO
693
  END IF
694

    
695

    
696

    
697
  !IF (IOptt .GE. 2) THEN
698
  !   Print  *, 'Counter=', Counter
699
  !END IF
700
  IF (Counter .GE. NbItMax) Then
701
     Print *, 'Counter GE 500, Stopping, normDeltaX_int=', normDeltaX_int
702
     STOP
703
  END IF
704

    
705

    
706
     x(:) = X_cart_k(1,:)
707
     y(:) = X_cart_k(2,:)
708
     z(:) = X_cart_k(3,:)
709

    
710

    
711
!  call CalcRmsd(Nat,x_k,y_k,z_k,x,y,z,MRot,rmsd,.TRUE.,.TRUE.)
712

    
713

    
714
  if (debug) WRITE(*,*) "COnverted in ",counter," iterations"
715

    
716
  DEALLOCATE(Geom,DeltaX_int,DeltaX_cart)
717
  DEALLOCATE(BprimT,BBT,BBT_inv,UMat_tmp)
718
  DEALLOCATE(X_cart_k)
719
  DEALLOCATE(GMat,EigVal,EigVec)
720

    
721
  XPrim=Xprimitive_t
722

    
723
  if (debug) WRITE(*,*) '=====ConvertBakerInternal_cart Over====='
724

    
725
END SUBROUTINE ConvertBakerInternal_cart