root / src / ConvGrad_Cart2Int.f90 @ 8
History  View  Annotate  Download (4.3 kB)
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 minimumnorm solution to a real linear least 
22 
! * squares problem: 
23 
! * minimize 2norm( b  A*x ) 
24 
! * using the singular value decomposition (SVD) of A. A is an MbyN 
25 
! * matrix which may be rankdeficient. 
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 = 1d8, eps = 1d12 
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 