root / src / Step_GEDIIS_All.f90 @ 8
Historique | Voir | Annoter | Télécharger (21,64 ko)
1 |
! Geom = input parameter vector (Geometry), Grad = input gradient vector, HEAT is Energy(Geom) |
---|---|
2 |
SUBROUTINE Step_GEDIIS_All(NGeomF,IGeom,Step,Geom,Grad,HEAT,Hess,NCoord,allocation_flag,Tangent) |
3 |
!SUBROUTINE Step_GEDIIS(Geom_new,Geom,Grad,HEAT,Hess,NCoord,FRST) |
4 |
use Io_module |
5 |
use Path_module, only : Nom, Atome, OrderInv, indzmat, Pi, Nat, Vfree |
6 |
IMPLICIT NONE |
7 |
|
8 |
INTEGER(KINT) :: NGeomF,IGeom |
9 |
INTEGER(KINT), INTENT(IN) :: NCoord |
10 |
REAL(KREAL) :: Geom(NCoord), Grad(NCoord), Hess(NCoord*NCoord), Step(NCoord) |
11 |
REAL(KREAL) :: HEAT ! HEAT= Energy |
12 |
LOGICAL :: allocation_flag |
13 |
REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord) |
14 |
|
15 |
! MRESET = maximum number of iterations. |
16 |
INTEGER(KINT), PARAMETER :: MRESET=15, M2=(MRESET+1)*(MRESET+1) !M2 = 256 |
17 |
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:,:), GradSet(:,:) ! NGeomF,MRESET*NCoord |
18 |
REAL(KREAL), ALLOCATABLE, SAVE :: GSAVE(:,:) !NGeomF,NCoord |
19 |
REAL(KREAL), ALLOCATABLE, SAVE :: ESET(:,:) |
20 |
REAL(KREAL) :: ESET_tmp(MRESET), B(M2), BS(M2), BST(M2), B_tmp(M2) ! M2=256 |
21 |
LOGICAL :: DEBUG, PRINT, ci_lt_zero |
22 |
INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:) ! mth Iteration |
23 |
LOGICAL, ALLOCATABLE, SAVE :: FRST(:) |
24 |
REAL(KREAL) :: ci(MRESET), ci_tmp(MRESET) ! MRESET = maximum number of iterations. |
25 |
INTEGER(KINT) :: NGEDIIS, MPLUS, INV, ITERA, MM, cis_zero |
26 |
INTEGER(KINT) :: I, J, K, JJ, JNV, II, IONE, IJ, IX, JX, KX |
27 |
INTEGER(KINT) :: current_size_B_mat, MyPointer, Isch, NFree, Idx |
28 |
REAL(KREAL) :: XMax, XNorm, DET, THRES, tmp, ER_star, ER_star_tmp, Norm |
29 |
REAL(KREAL), PARAMETER :: eps=1e-12 |
30 |
REAL(KREAL), PARAMETER :: crit=1e-8 |
31 |
REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord |
32 |
REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree |
33 |
REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree |
34 |
REAL(KREAL), ALLOCATABLE :: Grad_free(:), Step_free(:) ! NFree |
35 |
REAL(KREAL), ALLOCATABLE :: Geom_free(:), Geom_new_free(:) ! NFree |
36 |
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:), GradSet_free(:,:) |
37 |
|
38 |
DEBUG=.TRUE. |
39 |
PRINT=.FALSE. |
40 |
|
41 |
IF (PRINT) WRITE(*,'(/,'' BEGIN Step_GEDIIS_ALL '')') |
42 |
|
43 |
! Initialization |
44 |
IF (allocation_flag) THEN |
45 |
! allocation_flag will be set to False in SPACE_GEDIIS, so no need to modify it here |
46 |
IF (ALLOCATED(GeomSet)) THEN |
47 |
IF (PRINT) WRITE(*,'(/,'' In allocation_flag, GEDIIS_ALL Dealloc '')') |
48 |
DEALLOCATE(GeomSet,GradSet,GSave,GeomSet_free,GradSet_free) |
49 |
RETURN |
50 |
ELSE |
51 |
IF (PRINT) WRITE(*,'(/,'' In allocation_flag, GEDIIS_ALL Alloc '')') |
52 |
ALLOCATE(GeomSet(NGeomF,MRESET*NCoord),GradSet(NGeomF,MRESET*NCoord),GSAVE(NGeomF,NCoord)) |
53 |
ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord)) |
54 |
ALLOCATE(MSET(NGeomF),FRST(NGeomF),ESET(NGeomF,MRESET)) |
55 |
DO I=1,NGeomF |
56 |
FRST(I) = .TRUE. |
57 |
GeomSet(I,:) = 0.d0 |
58 |
GradSet(I,:) = 0.d0 |
59 |
GSAVE(I,:)=0.d0 |
60 |
GeomSet_free(I,:) = 0.d0 |
61 |
GradSet_free(I,:) = 0.d0 |
62 |
END DO |
63 |
MSET(:)=0 |
64 |
END IF |
65 |
allocation_flag = .FALSE. |
66 |
END IF ! IF (allocation_flag) THEN |
67 |
|
68 |
! ADDED FROM HERE: |
69 |
Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord), as of now, an Identity matrix. |
70 |
! we orthogonalize Vfree to the tangent vector of this geom only if Tangent/=0.d0 |
71 |
Norm=sqrt(dot_product(Tangent,Tangent)) |
72 |
IF (Norm.GT.eps) THEN |
73 |
ALLOCATE(Tanf(NCoord)) |
74 |
|
75 |
! We normalize Tangent |
76 |
Tangent=Tangent/Norm |
77 |
|
78 |
! We convert Tangent into Vfree only displacements. This is useless for now (2007.Apr.23) |
79 |
! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints. |
80 |
DO I=1,NCoord |
81 |
Tanf(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Tangent) |
82 |
END DO |
83 |
Tangent=0.d0 |
84 |
DO I=1,NCoord |
85 |
Tangent=Tangent+Tanf(I)*Vfree(:,I) |
86 |
END DO |
87 |
! first we subtract Tangent from vfree |
88 |
DO I=1,NCoord |
89 |
Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent) |
90 |
Vfree(:,I)=Vfree(:,I)-Norm*Tangent |
91 |
END DO |
92 |
|
93 |
Idx=0 |
94 |
! Schmidt orthogonalization of the Vfree vectors |
95 |
DO I=1,NCoord |
96 |
! We subtract the first vectors, we do it twice as the Schmidt procedure is not numerically stable. |
97 |
DO Isch=1,2 |
98 |
DO J=1,Idx |
99 |
Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,J),(/NCoord/))) |
100 |
Vfree(:,I)=Vfree(:,I)-Norm*Vfree(:,J) |
101 |
END DO |
102 |
END DO |
103 |
Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,I),(/NCoord/))) |
104 |
IF (Norm.GE.crit) THEN |
105 |
Idx=Idx+1 |
106 |
Vfree(:,Idx)=Vfree(:,I)/sqrt(Norm) |
107 |
END IF |
108 |
END DO |
109 |
|
110 |
IF (Idx/= NCoord-1) THEN |
111 |
WRITE(*,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom |
112 |
WRITE(IOOut,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom |
113 |
STOP |
114 |
END IF |
115 |
|
116 |
DEALLOCATE(Tanf) |
117 |
NFree=Idx |
118 |
ELSE ! Tangent =0, matches IF (Norm.GT.eps) THEN |
119 |
if (debug) WRITE(*,*) "Tangent=0, using full displacement" |
120 |
NFree=NCoord |
121 |
END IF !IF (Norm.GT.eps) THEN |
122 |
|
123 |
if (debug) WRITE(*,*) 'DBG Step_GEDIIS_All, IGeom, NFree=', IGeom, NFree |
124 |
|
125 |
! We now calculate the new step |
126 |
! we project the hessian onto the free vectors |
127 |
ALLOCATE(HFree(NFree*NFree),Htmp(NCoord,NFree),Grad_free(NFree)) |
128 |
ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free(NFree)) |
129 |
DO J=1,NFree |
130 |
DO I=1,NCoord |
131 |
Htmp(I,J)=0.d0 |
132 |
DO K=1,NCoord |
133 |
Htmp(I,J)=Htmp(I,J)+Hess(((I-1)*NCoord)+K)*Vfree(K,J) |
134 |
END DO |
135 |
END DO |
136 |
END DO |
137 |
DO J=1,NFree |
138 |
DO I=1,NFree |
139 |
HFree(I+((J-1)*NFree))=0.d0 |
140 |
DO K=1,NCoord |
141 |
HFree(I+((J-1)*NFree))=HFree(I+((J-1)*NFree))+Vfree(K,I)*Htmp(K,J) |
142 |
END DO |
143 |
END DO |
144 |
END DO |
145 |
|
146 |
DO I=1,NFree |
147 |
Grad_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Grad) |
148 |
Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom) |
149 |
END DO |
150 |
!ADDED ENDS HERE.*********************************************** |
151 |
|
152 |
! SPACE_GEDIIS SIMPLY LOADS THE CURRENT VALUES OF Geom AND Grad INTO THE ARRAYS GeomSet |
153 |
! AND GradSet, MSET is set to zero and then 1 in SPACE_GEDIIS_All at first iteration. |
154 |
CALL SPACE_GEDIIS_All(NGeomF,IGeom,MRESET,MSET,Geom,Grad,HEAT,NCoord,GeomSet,GradSet,ESET,FRST) |
155 |
|
156 |
IF (PRINT) WRITE(*,'(/,'' GEDIIS after SPACE_GEDIIS_ALL '')') |
157 |
|
158 |
DO J=1,MSet(IGeom) |
159 |
DO K=1,NFree |
160 |
GradSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),& |
161 |
GradSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord)) |
162 |
GeomSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),& |
163 |
GeomSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord)) |
164 |
END DO |
165 |
END DO |
166 |
|
167 |
! INITIALIZE SOME VARIABLES AND CONSTANTS: |
168 |
NGEDIIS = MSET(IGeom) !MSET=mth iteration |
169 |
MPLUS = MSET(IGeom) + 1 |
170 |
MM = MPLUS * MPLUS |
171 |
|
172 |
! CONSTRUCT THE GEDIIS MATRIX: |
173 |
! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)> |
174 |
JJ=0 |
175 |
INV=-NFree |
176 |
DO I=1,MSET(IGeom) |
177 |
INV=INV+NFree |
178 |
JNV=-NFree |
179 |
DO J=1,MSET(IGeom) |
180 |
JNV=JNV+NFree |
181 |
JJ = JJ + 1 |
182 |
B(JJ)=0.D0 |
183 |
DO K=1, NFree |
184 |
B(JJ) = B(JJ) + (((GradSet_free(IGeom,INV+K)-GradSet_free(IGeom,JNV+K))* & |
185 |
(GeomSet_free(IGeom,INV+K)-GeomSet_free(IGeom,JNV+K)))/2.D0) |
186 |
END DO |
187 |
END DO |
188 |
END DO |
189 |
|
190 |
! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix. |
191 |
! The correction is needed because the last coloumn of the matrix contains all 1 and one zero. |
192 |
DO I=MSET(IGeom)-1,1,-1 |
193 |
DO J=MSET(IGeom),1,-1 |
194 |
B(I*MSET(IGeom)+J+I) = B(I*MSET(IGeom)+J) |
195 |
END DO |
196 |
END DO |
197 |
|
198 |
! For the last row and last column of GEDIIS matrix: |
199 |
DO I=1,MPLUS |
200 |
B(MPLUS*I) = 1.D0 |
201 |
B(MPLUS*MSET(IGeom)+I) = 1.D0 |
202 |
END DO |
203 |
B(MM) = 0.D0 |
204 |
|
205 |
DO I=1, MPLUS |
206 |
!WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS)) |
207 |
END DO |
208 |
|
209 |
! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM: |
210 |
80 CONTINUE |
211 |
DO I=1,MM !MM = (MSET(IGeom)+1) * (MSET(IGeom)+1) |
212 |
BS(I) = B(I) !just a copy of the original B (GEDIIS) matrix |
213 |
END DO |
214 |
|
215 |
IF (NGEDIIS .NE. MSET(IGeom)) THEN |
216 |
DO II=1,MSET(IGeom)-NGEDIIS |
217 |
XMAX = -1.D10 |
218 |
ITERA = 0 |
219 |
DO I=1,MSET(IGeom) |
220 |
XNORM = 0.D0 |
221 |
INV = (I-1) * MPLUS |
222 |
DO J=1,MSET(IGeom) |
223 |
XNORM = XNORM + ABS(B(INV + J)) |
224 |
END DO |
225 |
IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN |
226 |
XMAX = XNORM |
227 |
ITERA = I |
228 |
IONE = INV + I |
229 |
ENDIF |
230 |
END DO |
231 |
|
232 |
DO I=1,MPLUS |
233 |
INV = (I-1) * MPLUS |
234 |
DO J=1,MPLUS |
235 |
JNV = (J-1) * MPLUS |
236 |
IF (J.EQ.ITERA) B(INV + J) = 0.D0 |
237 |
B(JNV + I) = B(INV + J) |
238 |
END DO |
239 |
END DO |
240 |
B(IONE) = 1.0D0 |
241 |
END DO |
242 |
END IF ! matches IF (NGEDIIS .NE. MSET(IGeom)) THEN |
243 |
|
244 |
! SCALE GEDIIS MATRIX BEFORE INVERSION: |
245 |
DO I=1,MPLUS |
246 |
II = MPLUS * (I-1) + I ! B(II)=diagonal elements of B matrix |
247 |
GSAVE(IGeom,I) = 1.D0 / DSQRT(1.D-20+DABS(B(II))) |
248 |
!Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I) |
249 |
END DO |
250 |
GSAVE(IGeom,MPLUS) = 1.D0 |
251 |
DO I=1,MPLUS |
252 |
DO J=1,MPLUS |
253 |
IJ = MPLUS * (I-1) + J |
254 |
B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J) |
255 |
END DO |
256 |
END DO |
257 |
|
258 |
! INVERT THE GEDIIS MATRIX B: |
259 |
DO I=1, MPLUS |
260 |
!WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS)) |
261 |
END DO |
262 |
|
263 |
CALL MINV(B,MPLUS,DET) ! matrix inversion. |
264 |
|
265 |
DO I=1, MPLUS |
266 |
!WRITE(*,'(10(1X,F20.16))') B((I-1)*MPLUS+1:I*(MPLUS)) |
267 |
END DO |
268 |
|
269 |
DO I=1,MPLUS |
270 |
DO J=1,MPLUS |
271 |
IJ = MPLUS * (I-1) + J |
272 |
B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J) |
273 |
END DO |
274 |
END DO |
275 |
|
276 |
! COMPUTE THE NEW INTERPOLATED PARAMETER VECTOR (Geometry): |
277 |
ci=0.d0 |
278 |
ci_tmp=0.d0 |
279 |
|
280 |
ci_lt_zero= .FALSE. |
281 |
DO I=1, MSET(IGeom) |
282 |
DO J=1, MSET(IGeom) ! B matrix is read column-wise |
283 |
ci(I)=ci(I)+B((J-1)*(MPLUS)+I)*ESET(IGeom,J) !ESET is energy set. |
284 |
END DO |
285 |
ci(I)=ci(I)+B((MPLUS-1)*(MPLUS)+I) |
286 |
!Print *, 'NO ci < 0 yet, c(',I,')=', ci(I) |
287 |
IF((ci(I) .LT. 0.0D0) .OR. (ci(I) .GT. 1.0D0)) THEN |
288 |
ci_lt_zero=.TRUE. |
289 |
EXIT |
290 |
END IF |
291 |
END DO !matches DO I=1, MSET(IGeom) |
292 |
|
293 |
IF (ci_lt_zero) Then |
294 |
cis_zero = 0 |
295 |
ER_star = 0.D0 |
296 |
ER_star_tmp = 1e32 |
297 |
|
298 |
! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)>, Full B matrix created first and then rows and columns are removed. |
299 |
JJ=0 |
300 |
INV=-NFree |
301 |
DO IX=1,MSET(IGeom) |
302 |
INV=INV+NFree |
303 |
JNV=-NFree |
304 |
DO JX=1,MSET(IGeom) |
305 |
JNV=JNV+NFree |
306 |
JJ = JJ + 1 |
307 |
BST(JJ)=0.D0 |
308 |
DO KX=1, NFree |
309 |
BST(JJ) = BST(JJ) + (((GradSet_free(IGeom,INV+KX)-GradSet_free(IGeom,JNV+KX))* & |
310 |
(GeomSet_free(IGeom,INV+KX)-GeomSet_free(IGeom,JNV+KX)))/2.D0) |
311 |
END DO |
312 |
END DO |
313 |
END DO |
314 |
|
315 |
DO I=1, (2**MSET(IGeom))-2 ! all (2**MSET(IGeom))-2 combinations of cis, except the one where all cis are .GT. 0 and .LT. 1 |
316 |
!Print *, 'Entering into DO I=1, (2**MSET(IGeom))-2 loop, MSET(IGeom)=', MSET(IGeom), ', I=', I |
317 |
ci(:)=1.D0 |
318 |
! find out which cis are zero in I: |
319 |
DO IX=1, MSET(IGeom) |
320 |
JJ=IAND(I, 2**(IX-1)) |
321 |
IF(JJ .EQ. 0) Then |
322 |
ci(IX)=0.D0 |
323 |
END IF |
324 |
END DO |
325 |
|
326 |
ci_lt_zero = .FALSE. |
327 |
! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)>, Full B matrix created first and then rows and columns are removed. |
328 |
DO IX=1, MSET(IGeom)*MSET(IGeom) |
329 |
B(IX) = BST(IX) !just a copy of the original B (GEDIIS) matrix |
330 |
END DO |
331 |
|
332 |
! Removal of KXth row and KXth column in order to accomodate cI to be zero: |
333 |
current_size_B_mat=MSET(IGeom) |
334 |
cis_zero = 0 |
335 |
! The bits of I (index of the upper loop 'DO I=1, (2**MSET(IGeom))-2'), gives which cis are zero. |
336 |
DO KX=1, MSET(IGeom) ! searching for each bit of I (index of the upper loop 'DO I=1, (2**MSET(IGeom))-2') |
337 |
IF (ci(KX) .EQ. 0.D0) Then !remove KXth row and KXth column |
338 |
cis_zero = cis_zero + 1 |
339 |
|
340 |
! First row removal: (B matrix is read column-wise) |
341 |
JJ=0 |
342 |
DO IX=1,current_size_B_mat ! columns reading |
343 |
DO JX=1,current_size_B_mat ! rows reading |
344 |
IF (JX .NE. KX) Then |
345 |
JJ = JJ + 1 |
346 |
B_tmp(JJ) = B((IX-1)*current_size_B_mat+JX) |
347 |
END IF |
348 |
END DO |
349 |
END DO |
350 |
|
351 |
DO IX=1,current_size_B_mat*(current_size_B_mat-1) |
352 |
B(IX) = B_tmp(IX) |
353 |
END DO |
354 |
|
355 |
! Now column removal: |
356 |
JJ=0 |
357 |
DO IX=1,KX-1 ! columns reading |
358 |
DO JX=1,current_size_B_mat-1 ! rows reading |
359 |
JJ = JJ + 1 |
360 |
B_tmp(JJ) = B(JJ) |
361 |
END DO |
362 |
END DO |
363 |
|
364 |
DO IX=KX+1,current_size_B_mat |
365 |
DO JX=1,current_size_B_mat-1 |
366 |
JJ = JJ + 1 |
367 |
B_tmp(JJ) = B(JJ+current_size_B_mat-1) |
368 |
END DO |
369 |
END DO |
370 |
|
371 |
DO IX=1,(current_size_B_mat-1)*(current_size_B_mat-1) |
372 |
B(IX) = B_tmp(IX) |
373 |
END DO |
374 |
current_size_B_mat = current_size_B_mat - 1 |
375 |
END IF ! matches IF (ci(KX) .EQ. 0.D0) Then !remove |
376 |
END DO ! matches DO KX=1, MSET(IGeom) |
377 |
|
378 |
! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix. |
379 |
! The correction is needed because the last coloumn and row of the matrix contains all 1 and one zero. |
380 |
DO IX=MSET(IGeom)-cis_zero-1,1,-1 |
381 |
DO JX=MSET(IGeom)-cis_zero,1,-1 |
382 |
B(IX*(MSET(IGeom)-cis_zero)+JX+IX) = B(IX*(MSET(IGeom)-cis_zero)+JX) |
383 |
END DO |
384 |
END DO |
385 |
|
386 |
! for last row and last column of GEDIIS matrix |
387 |
DO IX=1,MPLUS-cis_zero |
388 |
B((MPLUS-cis_zero)*IX) = 1.D0 |
389 |
B((MPLUS-cis_zero)*(MSET(IGeom)-cis_zero)+IX) = 1.D0 |
390 |
END DO |
391 |
B((MPLUS-cis_zero) * (MPLUS-cis_zero)) = 0.D0 |
392 |
|
393 |
DO IX=1, MPLUS |
394 |
!WRITE(*,'(10(1X,F20.4))') B((IX-1)*MPLUS+1:IX*(MPLUS)) |
395 |
END DO |
396 |
|
397 |
! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM: |
398 |
IF (NGEDIIS .NE. MSET(IGeom)) THEN |
399 |
JX = min(MSET(IGeom)-NGEDIIS,MSET(IGeom)-cis_zero-1) |
400 |
DO II=1,JX |
401 |
XMAX = -1.D10 |
402 |
ITERA = 0 |
403 |
DO IX=1,MSET(IGeom)-cis_zero |
404 |
XNORM = 0.D0 |
405 |
INV = (IX-1) * (MPLUS-cis_zero) |
406 |
DO J=1,MSET(IGeom)-cis_zero |
407 |
XNORM = XNORM + ABS(B(INV + J)) |
408 |
END DO |
409 |
IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN |
410 |
XMAX = XNORM |
411 |
ITERA = IX |
412 |
IONE = INV + IX |
413 |
ENDIF |
414 |
END DO |
415 |
|
416 |
DO IX=1,MPLUS-cis_zero |
417 |
INV = (IX-1) * (MPLUS-cis_zero) |
418 |
DO J=1,MPLUS-cis_zero |
419 |
JNV = (J-1) * (MPLUS-cis_zero) |
420 |
IF (J.EQ.ITERA) B(INV + J) = 0.D0 |
421 |
B(JNV + IX) = B(INV + J) |
422 |
END DO |
423 |
END DO |
424 |
B(IONE) = 1.0D0 |
425 |
END DO |
426 |
END IF ! matches IF (NGEDIIS .NE. MSET(IGeom)) THEN |
427 |
|
428 |
! SCALE GEDIIS MATRIX BEFORE INVERSION: |
429 |
DO IX=1,MPLUS-cis_zero |
430 |
II = (MPLUS-cis_zero) * (IX-1) + IX ! B(II)=diagonal elements of B matrix |
431 |
GSAVE(IGeom,IX) = 1.D0 / DSQRT(1.D-20+DABS(B(II))) |
432 |
END DO |
433 |
GSAVE(IGeom,MPLUS-cis_zero) = 1.D0 |
434 |
DO IX=1,MPLUS-cis_zero |
435 |
DO JX=1,MPLUS-cis_zero |
436 |
IJ = (MPLUS-cis_zero) * (IX-1) + JX |
437 |
B(IJ) = B(IJ) * GSAVE(IGeom,IX) * GSAVE(IGeom,JX) |
438 |
END DO |
439 |
END DO |
440 |
|
441 |
! INVERT THE GEDIIS MATRIX B: |
442 |
CALL MINV(B,MPLUS-cis_zero,DET) ! matrix inversion. |
443 |
|
444 |
DO IX=1,MPLUS-cis_zero |
445 |
DO JX=1,MPLUS-cis_zero |
446 |
IJ = (MPLUS-cis_zero) * (IX-1) + JX |
447 |
B(IJ) = B(IJ) * GSAVE(IGeom,IX) * GSAVE(IGeom,JX) |
448 |
END DO |
449 |
END DO |
450 |
|
451 |
DO IX=1, MPLUS |
452 |
!WRITE(*,'(10(1X,F20.4))') B((IX-1)*MPLUS+1:IX*(MPLUS)) |
453 |
END DO |
454 |
|
455 |
! ESET is rearranged to handle zero cis and stored in ESET_tmp: |
456 |
JJ=0 |
457 |
DO IX=1, MSET(IGeom) |
458 |
IF (ci(IX) .NE. 0) Then |
459 |
JJ=JJ+1 |
460 |
ESET_tmp(JJ) = ESET(IGeom,IX) |
461 |
END IF |
462 |
END DO |
463 |
|
464 |
! DETERMINATION OF nonzero cis: |
465 |
MyPointer=1 |
466 |
DO IX=1, MSET(IGeom)-cis_zero |
467 |
tmp = 0.D0 |
468 |
DO J=1, MSET(IGeom)-cis_zero ! B matrix is read column-wise |
469 |
tmp=tmp+B((J-1)*(MPLUS-cis_zero)+IX)*ESET_tmp(J) |
470 |
END DO |
471 |
tmp=tmp+B((MPLUS-cis_zero-1)*(MPLUS-cis_zero)+IX) |
472 |
IF((tmp .LT. 0.0D0) .OR. (tmp .GT. 1.0D0)) THEN |
473 |
ci_lt_zero=.TRUE. |
474 |
EXIT |
475 |
ELSE |
476 |
DO JX=MyPointer,MSET(IGeom) |
477 |
IF (ci(JX) .NE. 0) Then |
478 |
ci(JX) = tmp |
479 |
MyPointer=JX+1 |
480 |
EXIT |
481 |
END IF |
482 |
END DO |
483 |
END IF |
484 |
END DO !matches DO I=1, MSET(IGeom)-cis_zero |
485 |
!Print *, 'Local set of cis, first 10:, MSET(IGeom)=', MSET(IGeom), ', I of (2**MSET(IGeom))-2=', I |
486 |
!WRITE(*,'(10(1X,F20.4))') ci(1:MSET(IGeom)) |
487 |
!Print *, 'Local set of cis ends:****************************************' |
488 |
|
489 |
! new set of cis determined based on the lower energy (ER_star): |
490 |
IF (.NOT. ci_lt_zero) Then |
491 |
Call Energy_GEDIIS(MRESET,MSET(IGeom),ci,GeomSet_free(IGeom,:),GradSet_free(IGeom,:),ESET(IGeom,:),NFree,ER_star) |
492 |
IF (ER_star .LT. ER_star_tmp) Then |
493 |
ci_tmp=ci |
494 |
ER_star_tmp = ER_star |
495 |
END IF |
496 |
END IF ! matches IF (.NOT. ci_lt_zero) Then |
497 |
END DO !matches DO I=1, (2**K)-2 ! all (2**K)-2 combinations of cis, except the one where all cis are .GT. 0 and .LT. 1 |
498 |
ci = ci_tmp |
499 |
END IF! matches IF (ci_lt_zero) Then |
500 |
|
501 |
Print *, 'Final set of cis, first 10:***********************************' |
502 |
WRITE(*,'(10(1X,F20.4))') ci(1:MSET(IGeom)) |
503 |
Print *, 'Final set of cis ends:****************************************' |
504 |
Geom_new_free(:) = 0.D0 |
505 |
DO I=1, MSET(IGeom) |
506 |
Geom_new_free(:) = Geom_new_free(:) + (ci(I)*GeomSet_free(IGeom,(I-1)*NFree+1:I*NFree)) !MPLUS=MSET(IGeom)+1 |
507 |
! R_(N+1)=R*+DeltaR: |
508 |
DO J=1, NFree |
509 |
tmp=0.D0 |
510 |
DO K=1,NFree |
511 |
! this can be commented: |
512 |
!tmp=tmp+HFree((J-1)*NFree+K)*GradSet_free(IGeom,(I-1)*NFree+K) ! If Hinv=.False., then we need to invert Hess |
513 |
END DO |
514 |
Geom_new_free(J) = Geom_new_free(J) - (ci(I)*tmp) |
515 |
END DO |
516 |
END DO |
517 |
|
518 |
Step_free(:) = Geom_new_free(:) - Geom_free(:) |
519 |
|
520 |
XNORM = SQRT(DOT_PRODUCT(Step_free,Step_free)) |
521 |
IF (PRINT) THEN |
522 |
WRITE (6,'(/10X,''DEVIATION IN X '',F10.4,8X,''DETERMINANT '',G9.3)') XNORM, DET |
523 |
!WRITE(*,'(10X,''GEDIIS COEFFICIENTS'')') |
524 |
!WRITE(*,'(10X,5F12.5)') (B(MPLUS*MSET(IGeom)+I),I=1,MSET(IGeom)) |
525 |
ENDIF |
526 |
|
527 |
! THE FOLLOWING TOLERENCES FOR XNORM AND DET ARE SOMEWHAT ARBITRARY! |
528 |
THRES = MAX(10.D0**(-NFree), 1.D-25) |
529 |
IF (XNORM.GT.2.D0 .OR. DABS(DET) .LT. THRES) THEN |
530 |
IF (PRINT)THEN |
531 |
WRITE(*,*) "THE GEDIIS MATRIX IS ILL CONDITIONED" |
532 |
WRITE(*,*) " - PROBABLY, VECTORS ARE LINEARLY DEPENDENT - " |
533 |
WRITE(*,*) "THE GEDIIS STEP WILL BE REPEATED WITH A SMALLER SPACE" |
534 |
END IF |
535 |
DO K=1,MM |
536 |
B(K) = BS(K) ! why this is reverted? Because "IF (NGEDIIS .GT. 0) GO TO 80", see below |
537 |
END DO |
538 |
NGEDIIS = NGEDIIS - 1 |
539 |
IF (NGEDIIS .GT. 0) GO TO 80 |
540 |
IF (PRINT) WRITE(*,'(10X,''NEWTON-RAPHSON STEP TAKEN'')') |
541 |
Geom_new_free(:) = Geom_free(:) ! Geom_new is set to original Geom, thus Step = Geom(:) - Geom_new(:)=zero, the whole |
542 |
! new update to Geom_new is discarded, since XNORM.GT.2.D0 .OR. DABS(DET) .LT. THRES |
543 |
END IF ! matches IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN |
544 |
|
545 |
!****************************************************************************************************************** |
546 |
Geom_new_free(:) = 0.D0 |
547 |
DO I=1, MSET(IGeom) |
548 |
Geom_new_free(:) = Geom_new_free(:) + (ci(I)*GeomSet_free(IGeom,(I-1)*NFree+1:I*NFree)) !MPLUS=MSET(IGeom)+1 |
549 |
! R_(N+1)=R*+DeltaR: |
550 |
DO J=1, NFree |
551 |
tmp=0.D0 |
552 |
DO K=1,NFree |
553 |
tmp=tmp+HFree((J-1)*NFree+K)*GradSet_free(IGeom,(I-1)*NFree+K) ! If Hinv=.False., then we need to invert Hess |
554 |
END DO |
555 |
Geom_new_free(J) = Geom_new_free(J) - (ci(I)*tmp) |
556 |
END DO |
557 |
END DO |
558 |
|
559 |
Step_free(:) = Geom_new_free(:) - Geom_free(:) |
560 |
!****************************************************************************************************************** |
561 |
Step = 0.d0 |
562 |
DO I=1,NFree |
563 |
Step = Step + Step_free(I)*Vfree(:,I) |
564 |
END DO |
565 |
|
566 |
DEALLOCATE(Hfree,Htmp,Grad_free,Step_free,Geom_free,Geom_new_free) |
567 |
|
568 |
IF (PRINT) WRITE(*,'(/,'' END Step_GEDIIS_ALL '',/)') |
569 |
|
570 |
END SUBROUTINE Step_GEDIIS_All |
571 |
|