Statistiques
| Révision :

root / src / ConvGrad_Cart2Int.f90 @ 8

Historique | Voir | Annoter | Télécharger (4,29 ko)

1
      SUBROUTINE ConvGrad_Cart2Int(n,BMat,indzmat,GradCart,GradInt)
2

    
3
        use Io_module
4
        use Path_module, only : atome
5

    
6
        IMPLICIT NONE
7

    
8
!
9
        INTEGER(KINT), INTENT(IN) ::  n, indzmat(n,5)
10
        REAL(KREAL), INTENT(IN)  :: BMat(3,n,3,n), GradCart(3,n)
11
        REAL(KREAL), INTENT(INOUT)  :: GradInt(3,n)
12

    
13

    
14
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15
!
16
! This routines converts a gradient in cartesian coordinates into
17
! the gradient in internal coordinates (either Zmat, Mixed or Baker)
18
!
19
! It uses a LAPACK subroutine (dgelsd) which purpose is:
20

    
21
! *  DGELSD computes the minimum-norm solution to a real linear least
22
! *  squares problem:
23
! *      minimize 2-norm(| b - A*x |)
24
! *  using the singular value decomposition (SVD) of A. A is an M-by-N
25
! *  matrix which may be rank-deficient.
26

    
27
! Variables:
28
! input  : n, number of atoms
29
!          BMat: the dZ/dx matrix
30
!          GradCart: Gradient in cartesian coordinates
31
!
32
! output : GradInt: Gradient in internal coordintaes
33
!
34
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
35

    
36
      logical    debug, init
37
      REAL(KREAL), parameter :: zero = 0d0, tol = 1d-8,  eps = 1d-12
38
 
39
      INTEGER(KINT) :: n3, iat, k
40
      INTEGER(KINT) :: i
41
!
42
      save  debug, init
43
      data  init /.false./
44

    
45
!     .. Local Arrays  for LAPACK
46
      INTEGER(KINT) :: Info, LWork, Rank
47
      REAL(KREAL), ALLOCATABLE :: Work(:), GradCTmp(:, :)
48

    
49
!     IWORK dimension should be at least 3*MIN(M,N)*NLVL + 11*MIN(M,N),
50
!     where NLVL = MAX( 0, INT( LOG_2( MIN(M,N)/(SMLSIZ+1) ) )+1 )
51
!     and SMLSIZ = 25 
52
      INTEGER(KINT) :: LIWork, NLVL, SMLSIZ=12 ! We take smaller value just to be sure
53
      INTEGER(KINT), ALLOCATABLE :: IWork(:) ! IWork
54

    
55

    
56
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
57
!
58

    
59
  INTERFACE
60
     function valid(string) result (isValid)
61
       CHARACTER(*), intent(in) :: string
62
       logical                  :: isValid
63
     END function VALID
64

    
65
     SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, &
66
                        WORK, LWORK, IWORK, INFO )
67
! *
68
! *  -- LAPACK driver routine (version 3.2.2) --
69
! *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
70
! *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
71
! *     June 2010
72
! *
73
! *     .. Scalar Arguments ..
74
      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
75
      DOUBLE PRECISION   RCOND
76
!*     ..
77
!*     .. Array Arguments ..
78
      INTEGER            IWORK( * )
79
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
80

    
81
    END SUBROUTINE DGELSD
82

    
83

    
84

    
85
     
86
  END INTERFACE
87

    
88
!
89
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
90

    
91

    
92
      if (.not.init) then
93
        init = .true.
94
        debug = valid ('GEOSTEP').OR.valid('ConvGrad')
95
      endif
96

    
97
!
98
      Allocate(GradCTmp(3,n))
99

    
100
      if (debug) WRITE(*,*) "================ Entering ConvGrad_Cart2Int ============="
101
!
102
      n3 = 3 * n
103
!
104
      
105
      NLVL= MAX( 0, INT( LOG( n3*1.d0/(SMLSIZ+1))/LOG(2.d0) ) +1 )
106
      LIWork=Max(1,3*n3*NLVL + 11*n3+1)
107

    
108
      Allocate(IWork(LIWork),WOrk(1))
109

    
110
      LWork=-1
111
      GradInt=0.d0
112

    
113
        call DGELSD( n3, n3, 1, BMat, n3, GradCart, n3, GradInt, Tol, RANK, &
114
                        WORK, LWORK, IWORK, INFO )
115

    
116
        LWork=int(Work(1))
117
        if (debug) WRITE(*,*) "DEBUG CONVGRAD_CART2INT LWORK,N3=",LWORK,N3
118

    
119
        DEALLOCATE(Work)
120
        ALLOCATE(Work(LWork))
121

    
122
        call DGELSD( n3, n3, 1, BMat, n3, GradCart, n3, GradInt, Tol, RANK, &
123
                        WORK, LWORK, IWORK, INFO )
124

    
125
        GradInt=GradCart
126
      if (debug) then
127
         WRITE(*,*) "DEBUG CONVGRAD_CART2INT -- GradInt after DGELSD"
128
         WRITE(*,*) "INFO=",INFO
129
         DO I=1,N
130
            WRITE(*,'(1X,I5,3(1X,F13.6))') I, GradInt(1:3,I)
131
         END DO
132
      END IF
133

    
134

    
135
! In case we used a Zmat, we might have dummy atoms.
136
! We are extra cautious and put their gradient components to 0
137

    
138
      do iat = 1, n
139
        do k = 1, 3
140
          if (atome(indzmat(iat,1)) == 0) then
141
            GradInt(k,iat) = zero
142
          endif
143
        end do
144
      end do
145

    
146
      where (abs(GradInt)<=eps) GradInt=zero
147

    
148
      if (debug) write (*,*) "DBG ConvGrad_Cart2Int GradInt=",GradInt
149

    
150
      DEALLOCATE(IWork,Work,GradCTmp)
151
!
152
!
153
      if (debug) WRITE(*,*) "================  ConvGrad_Cart2Int Over ============="
154
      return
155
! ======================================================================
156

    
157
    end SUBROUTINE ConvGrad_Cart2Int