1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <time.h>
5 #include <string.h>
6 #include "openmx_common.h"
7 #include "mpi.h"
8
9
10 static void Simple_Mixing_H(int MD_iter, int SCF_iter, int SCF_iter0 );
11 static void Pulay_Mixing_H(int MD_iter, int SCF_iter, int SCF_iter0 );
12 static void Pulay_Mixing_H_MultiSecant(int MD_iter, int SCF_iter, int SCF_iter0 );
13 static void Pulay_Mixing_H_with_One_Shot_Hessian(int MD_iter, int SCF_iter, int SCF_iter0 );
14 static void Inverse(int n, double **a, double **ia);
15
16
Mixing_H(int MD_iter,int SCF_iter,int SCF_iter0)17 double Mixing_H( int MD_iter, int SCF_iter, int SCF_iter0 )
18 {
19 double time0;
20 double TStime,TEtime;
21 int numprocs,myid,ID;
22
23 /* MPI */
24 MPI_Comm_size(mpi_comm_level1,&numprocs);
25 MPI_Comm_rank(mpi_comm_level1,&myid);
26
27 MPI_Barrier(mpi_comm_level1);
28 dtime(&TStime);
29
30 /*******************************************************
31 Simple mixing
32 *******************************************************/
33
34 if ( SCF_iter<=(Pulay_SCF-1) ){
35 Simple_Mixing_H( MD_iter, SCF_iter, SCF_iter0 );
36 }
37
38 /*******************************************************
39 Pulay's method:
40 Residual Minimazation Method (RMM) using
41 Direct Inversion in the Iterative Subspace (DIIS)
42 *******************************************************/
43
44 else{
45
46 Pulay_Mixing_H( MD_iter, SCF_iter, SCF_iter0 );
47
48 /*
49 Pulay_Mixing_H_with_One_Shot_Hessian( MD_iter, SCF_iter, SCF_iter0 );
50 */
51
52 /*
53 Pulay_Mixing_H_MultiSecant( MD_iter, SCF_iter, SCF_iter0 );
54 */
55 }
56
57 /* if SCF_iter0==1, then NormRD[0]=1 */
58 if (SCF_iter0==1) NormRD[0] = 1.0;
59
60 MPI_Barrier(mpi_comm_level1);
61 dtime(&TEtime);
62 time0 = TEtime - TStime;
63 return time0;
64 }
65
66
67
68
69
Pulay_Mixing_H_MultiSecant(int MD_iter,int SCF_iter,int SCF_iter0)70 void Pulay_Mixing_H_MultiSecant(int MD_iter, int SCF_iter, int SCF_iter0 )
71 {
72 int Mc_AN,Gc_AN,Cwan,Hwan,h_AN,Gh_AN,i,j,spin;
73 int dim,m,n,flag_nan;
74 double sum,my_sum,tmp1,tmp2,alpha;
75 double r,r10,r11,r12,r13,r20,r21,r22;
76 double h,h10,h11,h12,h13,h20,h21,h22;
77 double my_sy,my_yy,sy,yy,norm,s,y,orx,al,be;
78 double **A,**IA,*coes,*coes2,*ror;
79 char nanchar[300];
80
81 /****************************************************
82 determination of dimension of the subspace
83 ****************************************************/
84
85 if (SCF_iter<=Num_Mixing_pDM) dim = SCF_iter-1;
86 else dim = Num_Mixing_pDM;
87
88 /****************************************************
89 allocation of arrays
90 ****************************************************/
91
92 coes = (double*)malloc(sizeof(double)*List_YOUSO[39]);
93 coes2 = (double*)malloc(sizeof(double)*List_YOUSO[39]);
94 ror = (double*)malloc(sizeof(double)*List_YOUSO[39]);
95
96 A = (double**)malloc(sizeof(double*)*List_YOUSO[39]);
97 for (i=0; i<List_YOUSO[39]; i++){
98 A[i] = (double*)malloc(sizeof(double)*List_YOUSO[39]);
99 }
100
101 IA = (double**)malloc(sizeof(double*)*List_YOUSO[39]);
102 for (i=0; i<List_YOUSO[39]; i++){
103 IA[i] = (double*)malloc(sizeof(double)*List_YOUSO[39]);
104 }
105
106 /****************************************************
107 shift the residual H
108 ****************************************************/
109
110 if (SpinP_switch==3){
111
112 /* shift the residual Hamiltonian */
113
114 for (m=dim; 0<m; m--){
115 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
116 Gc_AN = M2G[Mc_AN];
117 Cwan = WhatSpecies[Gc_AN];
118 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
119 Gh_AN = natn[Gc_AN][h_AN];
120 Hwan = WhatSpecies[Gh_AN];
121 for (i=0; i<Spe_Total_NO[Cwan]; i++){
122 for (j=0; j<Spe_Total_NO[Hwan]; j++){
123
124 ResidualH1[m][0][Mc_AN][h_AN][i][j] = ResidualH1[m-1][0][Mc_AN][h_AN][i][j];
125 ResidualH1[m][1][Mc_AN][h_AN][i][j] = ResidualH1[m-1][1][Mc_AN][h_AN][i][j];
126 ResidualH1[m][2][Mc_AN][h_AN][i][j] = ResidualH1[m-1][2][Mc_AN][h_AN][i][j];
127 ResidualH1[m][3][Mc_AN][h_AN][i][j] = ResidualH1[m-1][3][Mc_AN][h_AN][i][j];
128
129 ResidualH2[m][0][Mc_AN][h_AN][i][j] = ResidualH2[m-1][0][Mc_AN][h_AN][i][j];
130 ResidualH2[m][1][Mc_AN][h_AN][i][j] = ResidualH2[m-1][1][Mc_AN][h_AN][i][j];
131 ResidualH2[m][2][Mc_AN][h_AN][i][j] = ResidualH2[m-1][2][Mc_AN][h_AN][i][j];
132 }
133 }
134 }
135 }
136 }
137
138 /* calculate the current residual Hamiltonian */
139
140 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
141 Gc_AN = M2G[Mc_AN];
142 Cwan = WhatSpecies[Gc_AN];
143 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
144 Gh_AN = natn[Gc_AN][h_AN];
145 Hwan = WhatSpecies[Gh_AN];
146 for (i=0; i<Spe_Total_NO[Cwan]; i++){
147 for (j=0; j<Spe_Total_NO[Hwan]; j++){
148
149 ResidualH1[0][0][Mc_AN][h_AN][i][j] = H[0][Mc_AN][h_AN][i][j] - HisH1[0][0][Mc_AN][h_AN][i][j];
150 ResidualH1[0][1][Mc_AN][h_AN][i][j] = H[1][Mc_AN][h_AN][i][j] - HisH1[0][1][Mc_AN][h_AN][i][j];
151 ResidualH1[0][2][Mc_AN][h_AN][i][j] = H[2][Mc_AN][h_AN][i][j] - HisH1[0][2][Mc_AN][h_AN][i][j];
152 ResidualH1[0][3][Mc_AN][h_AN][i][j] = H[3][Mc_AN][h_AN][i][j] - HisH1[0][3][Mc_AN][h_AN][i][j];
153
154 ResidualH2[0][0][Mc_AN][h_AN][i][j] = iHNL[0][Mc_AN][h_AN][i][j] - HisH2[0][0][Mc_AN][h_AN][i][j];
155 ResidualH2[0][1][Mc_AN][h_AN][i][j] = iHNL[1][Mc_AN][h_AN][i][j] - HisH2[0][1][Mc_AN][h_AN][i][j];
156 ResidualH2[0][2][Mc_AN][h_AN][i][j] = iHNL[2][Mc_AN][h_AN][i][j] - HisH2[0][2][Mc_AN][h_AN][i][j];
157 }
158 }
159 }
160 }
161
162 }
163
164 else{
165
166 /* shift the residual Hamiltonian */
167
168 for (m=dim; 0<m; m--){
169 for (spin=0; spin<=SpinP_switch; spin++){
170 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
171 Gc_AN = M2G[Mc_AN];
172 Cwan = WhatSpecies[Gc_AN];
173 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
174 Gh_AN = natn[Gc_AN][h_AN];
175 Hwan = WhatSpecies[Gh_AN];
176 for (i=0; i<Spe_Total_NO[Cwan]; i++){
177 for (j=0; j<Spe_Total_NO[Hwan]; j++){
178 ResidualH1[m][spin][Mc_AN][h_AN][i][j] = ResidualH1[m-1][spin][Mc_AN][h_AN][i][j];
179 }
180 }
181 }
182 }
183 }
184 }
185
186 /* calculate the current residual Hamiltonian */
187
188 for (spin=0; spin<=SpinP_switch; spin++){
189 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
190 Gc_AN = M2G[Mc_AN];
191 Cwan = WhatSpecies[Gc_AN];
192 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
193 Gh_AN = natn[Gc_AN][h_AN];
194 Hwan = WhatSpecies[Gh_AN];
195 for (i=0; i<Spe_Total_NO[Cwan]; i++){
196 for (j=0; j<Spe_Total_NO[Hwan]; j++){
197 ResidualH1[0][spin][Mc_AN][h_AN][i][j] = H[spin][Mc_AN][h_AN][i][j] - HisH1[0][spin][Mc_AN][h_AN][i][j];
198 }
199 }
200 }
201 }
202 }
203
204 } /* else */
205
206 /****************************************************
207 calculation of the residual matrix
208 ****************************************************/
209
210 for (m=0; m<dim; m++){
211 for (n=0; n<dim; n++){
212
213 my_sum = 0.0;
214
215 if (SpinP_switch==3){
216
217 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
218 Gc_AN = M2G[Mc_AN];
219 Cwan = WhatSpecies[Gc_AN];
220 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
221 Gh_AN = natn[Gc_AN][h_AN];
222 Hwan = WhatSpecies[Gh_AN];
223 for (i=0; i<Spe_Total_NO[Cwan]; i++){
224 for (j=0; j<Spe_Total_NO[Hwan]; j++){
225
226 tmp1 = ResidualH1[m][0][Mc_AN][h_AN][i][j];
227 tmp2 = ResidualH1[n][0][Mc_AN][h_AN][i][j];
228 my_sum += tmp1*tmp2;
229
230 tmp1 = ResidualH1[m][1][Mc_AN][h_AN][i][j];
231 tmp2 = ResidualH1[n][1][Mc_AN][h_AN][i][j];
232 my_sum += tmp1*tmp2;
233
234 tmp1 = ResidualH1[m][2][Mc_AN][h_AN][i][j];
235 tmp2 = ResidualH1[n][2][Mc_AN][h_AN][i][j];
236 my_sum += tmp1*tmp2;
237
238 tmp1 = ResidualH1[m][3][Mc_AN][h_AN][i][j];
239 tmp2 = ResidualH1[n][3][Mc_AN][h_AN][i][j];
240 my_sum += tmp1*tmp2;
241
242 tmp1 = ResidualH2[m][0][Mc_AN][h_AN][i][j];
243 tmp2 = ResidualH2[n][0][Mc_AN][h_AN][i][j];
244 my_sum += tmp1*tmp2;
245
246 tmp1 = ResidualH2[m][1][Mc_AN][h_AN][i][j];
247 tmp2 = ResidualH2[n][1][Mc_AN][h_AN][i][j];
248 my_sum += tmp1*tmp2;
249
250 tmp1 = ResidualH2[m][2][Mc_AN][h_AN][i][j];
251 tmp2 = ResidualH2[n][2][Mc_AN][h_AN][i][j];
252 my_sum += tmp1*tmp2;
253 }
254 }
255 }
256 }
257
258 } /* if (SpinP_switch==3 */
259
260 else{
261
262 for (spin=0; spin<=SpinP_switch; spin++){
263 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
264 Gc_AN = M2G[Mc_AN];
265 Cwan = WhatSpecies[Gc_AN];
266 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
267 Gh_AN = natn[Gc_AN][h_AN];
268 Hwan = WhatSpecies[Gh_AN];
269 for (i=0; i<Spe_Total_NO[Cwan]; i++){
270 for (j=0; j<Spe_Total_NO[Hwan]; j++){
271 tmp1 = ResidualH1[m][spin][Mc_AN][h_AN][i][j];
272 tmp2 = ResidualH1[n][spin][Mc_AN][h_AN][i][j];
273 my_sum += tmp1*tmp2;
274 }
275 }
276 }
277 }
278 }
279
280 } /* else */
281
282 MPI_Allreduce(&my_sum, &A[m][n], 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
283 A[n][m] = A[m][n];
284
285 } /* n */
286 } /* m */
287
288 NormRD[0] = A[0][0];
289
290 for (m=1; m<=dim; m++){
291 A[m-1][dim] = -1.0;
292 A[dim][m-1] = -1.0;
293 }
294 A[dim][dim] = 0.0;
295
296 Inverse(dim,A,IA);
297
298 for (m=1; m<=dim; m++){
299 coes[m] = -IA[m-1][dim];
300 }
301
302 /****************************************************
303 check "nan", "NaN", "inf" or "Inf"
304 ****************************************************/
305
306 flag_nan = 0;
307 for (m=1; m<=dim; m++){
308
309 sprintf(nanchar,"%8.4f",coes[m]);
310 if ( strstr(nanchar,"nan")!=NULL || strstr(nanchar,"NaN")!=NULL
311 || strstr(nanchar,"inf")!=NULL || strstr(nanchar,"Inf")!=NULL){
312
313 flag_nan = 1;
314 }
315 }
316
317 if (flag_nan==1){
318 for (m=1; m<=dim; m++){
319 coes[m] = 0.0;
320 }
321 coes[1] = 0.05;
322 coes[2] = 0.95;
323 }
324
325 /****************************************************
326 calculation of optimum residual Hamiltonian
327 ****************************************************/
328
329 if (SpinP_switch==3){
330
331 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
332 Gc_AN = M2G[Mc_AN];
333 Cwan = WhatSpecies[Gc_AN];
334 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
335 Gh_AN = natn[Gc_AN][h_AN];
336 Hwan = WhatSpecies[Gh_AN];
337 for (i=0; i<Spe_Total_NO[Cwan]; i++){
338 for (j=0; j<Spe_Total_NO[Hwan]; j++){
339
340 r10 = 0.0;
341 r11 = 0.0;
342 r12 = 0.0;
343 r13 = 0.0;
344
345 r20 = 0.0;
346 r21 = 0.0;
347 r22 = 0.0;
348
349 for (m=0; m<dim; m++){
350
351 r10 += ResidualH1[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
352 r11 += ResidualH1[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
353 r12 += ResidualH1[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
354 r13 += ResidualH1[m][3][Mc_AN][h_AN][i][j]*coes[m+1];
355
356 r20 += ResidualH2[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
357 r21 += ResidualH2[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
358 r22 += ResidualH2[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
359 }
360
361 /* optimum Residual H is stored in ResidualH1[dim+1] and ResidualH2[dim+1] */
362
363 ResidualH1[dim+1][0][Mc_AN][h_AN][i][j] = r10;
364 ResidualH1[dim+1][1][Mc_AN][h_AN][i][j] = r11;
365 ResidualH1[dim+1][2][Mc_AN][h_AN][i][j] = r12;
366 ResidualH1[dim+1][3][Mc_AN][h_AN][i][j] = r13;
367
368 ResidualH2[dim+1][0][Mc_AN][h_AN][i][j] = r20;
369 ResidualH2[dim+1][1][Mc_AN][h_AN][i][j] = r21;
370 ResidualH2[dim+1][2][Mc_AN][h_AN][i][j] = r22;
371
372 }
373 }
374 }
375 }
376
377 }
378
379 else{
380
381 for (spin=0; spin<=SpinP_switch; spin++){
382 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
383 Gc_AN = M2G[Mc_AN];
384 Cwan = WhatSpecies[Gc_AN];
385 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
386 Gh_AN = natn[Gc_AN][h_AN];
387 Hwan = WhatSpecies[Gh_AN];
388 for (i=0; i<Spe_Total_NO[Cwan]; i++){
389 for (j=0; j<Spe_Total_NO[Hwan]; j++){
390
391 r = 0.0;
392 for (m=0; m<dim; m++){
393 r += ResidualH1[m][spin][Mc_AN][h_AN][i][j]*coes[m+1];
394 }
395
396 /* optimum Residual H is stored in ResidualH1[dim+1] */
397
398 ResidualH1[dim+1][spin][Mc_AN][h_AN][i][j] = r;
399
400 }
401 }
402 }
403 }
404 }
405 }
406
407 /******************************************************
408 calculations of inner products of <s0|y0> and <y0|y0>
409 in order to estimate the parameter "al".
410 ******************************************************/
411
412 my_sy = 0.0;
413 my_yy = 0.0;
414
415 if (SpinP_switch==3){
416
417 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
418 Gc_AN = M2G[Mc_AN];
419 Cwan = WhatSpecies[Gc_AN];
420 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
421 Gh_AN = natn[Gc_AN][h_AN];
422 Hwan = WhatSpecies[Gh_AN];
423 for (i=0; i<Spe_Total_NO[Cwan]; i++){
424 for (j=0; j<Spe_Total_NO[Hwan]; j++){
425
426 tmp1 = HisH1[0][0][Mc_AN][h_AN][i][j] - HisH1[1][0][Mc_AN][h_AN][i][j]; /* s */
427 tmp2 = ResidualH1[0][0][Mc_AN][h_AN][i][j] - ResidualH1[1][0][Mc_AN][h_AN][i][j]; /* y */
428 my_sy += tmp1*tmp2;
429 my_yy += tmp2*tmp2;
430
431 tmp1 = HisH1[0][1][Mc_AN][h_AN][i][j] - HisH1[1][1][Mc_AN][h_AN][i][j]; /* s */
432 tmp2 = ResidualH1[0][1][Mc_AN][h_AN][i][j] - ResidualH1[1][1][Mc_AN][h_AN][i][j]; /* y */
433 my_sy += tmp1*tmp2;
434 my_yy += tmp2*tmp2;
435
436 tmp1 = HisH1[0][2][Mc_AN][h_AN][i][j] - HisH1[1][2][Mc_AN][h_AN][i][j]; /* s */
437 tmp2 = ResidualH1[0][2][Mc_AN][h_AN][i][j] - ResidualH1[1][2][Mc_AN][h_AN][i][j]; /* y */
438 my_sy += tmp1*tmp2;
439 my_yy += tmp2*tmp2;
440
441 tmp1 = HisH1[0][3][Mc_AN][h_AN][i][j] - HisH1[1][3][Mc_AN][h_AN][i][j]; /* s */
442 tmp2 = ResidualH1[0][3][Mc_AN][h_AN][i][j] - ResidualH1[1][3][Mc_AN][h_AN][i][j]; /* y */
443 my_sy += tmp1*tmp2;
444 my_yy += tmp2*tmp2;
445
446 tmp1 = HisH2[0][0][Mc_AN][h_AN][i][j] - HisH2[1][0][Mc_AN][h_AN][i][j]; /* s */
447 tmp2 = ResidualH2[0][0][Mc_AN][h_AN][i][j] - ResidualH2[1][0][Mc_AN][h_AN][i][j]; /* y */
448 my_sy += tmp1*tmp2;
449 my_yy += tmp2*tmp2;
450
451 tmp1 = HisH2[0][1][Mc_AN][h_AN][i][j] - HisH2[1][1][Mc_AN][h_AN][i][j]; /* s */
452 tmp2 = ResidualH2[0][1][Mc_AN][h_AN][i][j] - ResidualH2[1][1][Mc_AN][h_AN][i][j]; /* y */
453 my_sy += tmp1*tmp2;
454 my_yy += tmp2*tmp2;
455
456 tmp1 = HisH2[0][2][Mc_AN][h_AN][i][j] - HisH2[1][2][Mc_AN][h_AN][i][j]; /* s */
457 tmp2 = ResidualH2[0][2][Mc_AN][h_AN][i][j] - ResidualH2[1][2][Mc_AN][h_AN][i][j]; /* y */
458 my_sy += tmp1*tmp2;
459 my_yy += tmp2*tmp2;
460 }
461 }
462 }
463 }
464
465 } /* if (SpinP_switch==3 */
466
467 else{
468
469 for (spin=0; spin<=SpinP_switch; spin++){
470 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
471 Gc_AN = M2G[Mc_AN];
472 Cwan = WhatSpecies[Gc_AN];
473 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
474 Gh_AN = natn[Gc_AN][h_AN];
475 Hwan = WhatSpecies[Gh_AN];
476 for (i=0; i<Spe_Total_NO[Cwan]; i++){
477 for (j=0; j<Spe_Total_NO[Hwan]; j++){
478
479 tmp1 = HisH1[0][spin][Mc_AN][h_AN][i][j] - HisH1[1][spin][Mc_AN][h_AN][i][j]; /* s */
480 tmp2 = ResidualH1[0][spin][Mc_AN][h_AN][i][j] - ResidualH1[1][spin][Mc_AN][h_AN][i][j]; /* y */
481 my_sy += tmp1*tmp2;
482 my_yy += tmp2*tmp2;
483 }
484 }
485 }
486 }
487 }
488
489 } /* else */
490
491 MPI_Allreduce(&my_sy, &sy, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
492 MPI_Allreduce(&my_yy, &yy, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
493
494 /* al < sy/yy */
495
496 al = sy/yy - 0.2;
497
498 /****************************************************
499 calculations of inner products of <r|y>
500 ****************************************************/
501
502 for (m=0; m<dim; m++){
503 for (n=0; n<dim; n++){
504
505 my_sum = 0.0;
506
507 if (SpinP_switch==3){
508
509 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
510 Gc_AN = M2G[Mc_AN];
511 Cwan = WhatSpecies[Gc_AN];
512 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
513 Gh_AN = natn[Gc_AN][h_AN];
514 Hwan = WhatSpecies[Gh_AN];
515 for (i=0; i<Spe_Total_NO[Cwan]; i++){
516 for (j=0; j<Spe_Total_NO[Hwan]; j++){
517
518 /* m */
519 s = HisH1[m][0][Mc_AN][h_AN][i][j] - HisH1[m+1][0][Mc_AN][h_AN][i][j]; /* s */
520 y = ResidualH1[m][0][Mc_AN][h_AN][i][j] - ResidualH1[m+1][0][Mc_AN][h_AN][i][j]; /* y */
521 r = s - al*y; /* r */
522 /* n */
523 y = ResidualH1[n][0][Mc_AN][h_AN][i][j] - ResidualH1[n+1][0][Mc_AN][h_AN][i][j]; /* y */
524 my_sum += r*y;
525
526 /* m */
527 s = HisH1[m][1][Mc_AN][h_AN][i][j] - HisH1[m+1][1][Mc_AN][h_AN][i][j]; /* s */
528 y = ResidualH1[m][1][Mc_AN][h_AN][i][j] - ResidualH1[m+1][1][Mc_AN][h_AN][i][j]; /* y */
529 r = s - al*y; /* r */
530 /* n */
531 y = ResidualH1[n][1][Mc_AN][h_AN][i][j] - ResidualH1[n+1][1][Mc_AN][h_AN][i][j]; /* y */
532 my_sum += r*y;
533
534 /* m */
535 s = HisH1[m][2][Mc_AN][h_AN][i][j] - HisH1[m+1][2][Mc_AN][h_AN][i][j]; /* s */
536 y = ResidualH1[m][2][Mc_AN][h_AN][i][j] - ResidualH1[m+1][2][Mc_AN][h_AN][i][j]; /* y */
537 r = s - al*y; /* r */
538 /* n */
539 y = ResidualH1[n][2][Mc_AN][h_AN][i][j] - ResidualH1[n+1][2][Mc_AN][h_AN][i][j]; /* y */
540 my_sum += r*y;
541
542 /* m */
543 s = HisH1[m][3][Mc_AN][h_AN][i][j] - HisH1[m+1][3][Mc_AN][h_AN][i][j]; /* s */
544 y = ResidualH1[m][3][Mc_AN][h_AN][i][j] - ResidualH1[m+1][3][Mc_AN][h_AN][i][j]; /* y */
545 r = s - al*y; /* r */
546 /* n */
547 y = ResidualH1[n][3][Mc_AN][h_AN][i][j] - ResidualH1[n+1][3][Mc_AN][h_AN][i][j]; /* y */
548 my_sum += r*y;
549
550 /* m */
551 s = HisH2[m][0][Mc_AN][h_AN][i][j] - HisH2[m+1][0][Mc_AN][h_AN][i][j]; /* s */
552 y = ResidualH2[m][0][Mc_AN][h_AN][i][j] - ResidualH2[m+1][0][Mc_AN][h_AN][i][j]; /* y */
553 r = s - al*y; /* r */
554 /* n */
555 y = ResidualH2[n][0][Mc_AN][h_AN][i][j] - ResidualH2[n+1][0][Mc_AN][h_AN][i][j]; /* y */
556 my_sum += r*y;
557
558 /* m */
559 s = HisH2[m][1][Mc_AN][h_AN][i][j] - HisH2[m+1][1][Mc_AN][h_AN][i][j]; /* s */
560 y = ResidualH2[m][1][Mc_AN][h_AN][i][j] - ResidualH2[m+1][1][Mc_AN][h_AN][i][j]; /* y */
561 r = s - al*y; /* r */
562 /* n */
563 y = ResidualH2[n][1][Mc_AN][h_AN][i][j] - ResidualH2[n+1][1][Mc_AN][h_AN][i][j]; /* y */
564 my_sum += r*y;
565
566 /* m */
567 s = HisH2[m][2][Mc_AN][h_AN][i][j] - HisH2[m+1][2][Mc_AN][h_AN][i][j]; /* s */
568 y = ResidualH2[m][2][Mc_AN][h_AN][i][j] - ResidualH2[m+1][2][Mc_AN][h_AN][i][j]; /* y */
569 r = s - al*y; /* r */
570 /* n */
571 y = ResidualH2[n][2][Mc_AN][h_AN][i][j] - ResidualH2[n+1][2][Mc_AN][h_AN][i][j]; /* y */
572 my_sum += r*y;
573
574 }
575 }
576 }
577 }
578
579 } /* if (SpinP_switch==3 */
580
581 else{
582
583 for (spin=0; spin<=SpinP_switch; spin++){
584 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
585 Gc_AN = M2G[Mc_AN];
586 Cwan = WhatSpecies[Gc_AN];
587 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
588 Gh_AN = natn[Gc_AN][h_AN];
589 Hwan = WhatSpecies[Gh_AN];
590 for (i=0; i<Spe_Total_NO[Cwan]; i++){
591 for (j=0; j<Spe_Total_NO[Hwan]; j++){
592
593 /* m */
594 s = HisH1[m][spin][Mc_AN][h_AN][i][j] - HisH1[m+1][spin][Mc_AN][h_AN][i][j]; /* s */
595 y = ResidualH1[m][spin][Mc_AN][h_AN][i][j] - ResidualH1[m+1][spin][Mc_AN][h_AN][i][j]; /* y */
596 r = s - al*y; /* r */
597 /* n */
598 y = ResidualH1[n][spin][Mc_AN][h_AN][i][j] - ResidualH1[n+1][spin][Mc_AN][h_AN][i][j]; /* y */
599 my_sum += r*y;
600
601 }
602 }
603 }
604 }
605 }
606
607 } /* else */
608
609 MPI_Allreduce(&my_sum, &A[m][n], 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
610
611 } /* n */
612 } /* m */
613
614 Inverse(dim-1,A,IA);
615
616 /****************************************************
617 calculations of inner products of <r|OptResidualH>
618 ****************************************************/
619
620 for (m=0; m<dim; m++){
621
622 my_sum = 0.0;
623
624 if (SpinP_switch==3){
625
626 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
627 Gc_AN = M2G[Mc_AN];
628 Cwan = WhatSpecies[Gc_AN];
629 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
630 Gh_AN = natn[Gc_AN][h_AN];
631 Hwan = WhatSpecies[Gh_AN];
632 for (i=0; i<Spe_Total_NO[Cwan]; i++){
633 for (j=0; j<Spe_Total_NO[Hwan]; j++){
634
635 s = HisH1[m][0][Mc_AN][h_AN][i][j] - HisH1[m+1][0][Mc_AN][h_AN][i][j]; /* s */
636 y = ResidualH1[m][0][Mc_AN][h_AN][i][j] - ResidualH1[m+1][0][Mc_AN][h_AN][i][j]; /* y */
637 r = s - al*y; /* r */
638 orx = ResidualH1[dim+1][0][Mc_AN][h_AN][i][j]; /* OptResidualH */
639 my_sum += r*orx;
640
641 s = HisH1[m][1][Mc_AN][h_AN][i][j] - HisH1[m+1][1][Mc_AN][h_AN][i][j]; /* s */
642 y = ResidualH1[m][1][Mc_AN][h_AN][i][j] - ResidualH1[m+1][1][Mc_AN][h_AN][i][j]; /* y */
643 r = s - al*y; /* r */
644 orx = ResidualH1[dim+1][1][Mc_AN][h_AN][i][j]; /* OptResidualH */
645 my_sum += r*orx;
646
647 s = HisH1[m][2][Mc_AN][h_AN][i][j] - HisH1[m+1][2][Mc_AN][h_AN][i][j]; /* s */
648 y = ResidualH1[m][2][Mc_AN][h_AN][i][j] - ResidualH1[m+1][2][Mc_AN][h_AN][i][j]; /* y */
649 r = s - al*y; /* r */
650 orx = ResidualH1[dim+1][2][Mc_AN][h_AN][i][j]; /* OptResidualH */
651 my_sum += r*orx;
652
653 s = HisH1[m][3][Mc_AN][h_AN][i][j] - HisH1[m+1][3][Mc_AN][h_AN][i][j]; /* s */
654 y = ResidualH1[m][3][Mc_AN][h_AN][i][j] - ResidualH1[m+1][3][Mc_AN][h_AN][i][j]; /* y */
655 r = s - al*y; /* r */
656 orx = ResidualH1[dim+1][3][Mc_AN][h_AN][i][j]; /* OptResidualH */
657 my_sum += r*orx;
658
659 s = HisH2[m][0][Mc_AN][h_AN][i][j] - HisH2[m+1][0][Mc_AN][h_AN][i][j]; /* s */
660 y = ResidualH2[m][0][Mc_AN][h_AN][i][j] - ResidualH2[m+1][0][Mc_AN][h_AN][i][j]; /* y */
661 r = s - al*y; /* r */
662 orx = ResidualH2[dim+1][0][Mc_AN][h_AN][i][j]; /* OptResidualH */
663 my_sum += r*orx;
664
665 s = HisH2[m][1][Mc_AN][h_AN][i][j] - HisH2[m+1][1][Mc_AN][h_AN][i][j]; /* s */
666 y = ResidualH2[m][1][Mc_AN][h_AN][i][j] - ResidualH2[m+1][1][Mc_AN][h_AN][i][j]; /* y */
667 r = s - al*y; /* r */
668 orx = ResidualH2[dim+1][1][Mc_AN][h_AN][i][j]; /* OptResidualH */
669 my_sum += r*orx;
670
671 s = HisH2[m][2][Mc_AN][h_AN][i][j] - HisH2[m+1][2][Mc_AN][h_AN][i][j]; /* s */
672 y = ResidualH2[m][2][Mc_AN][h_AN][i][j] - ResidualH2[m+1][2][Mc_AN][h_AN][i][j]; /* y */
673 r = s - al*y; /* r */
674 orx = ResidualH2[dim+1][2][Mc_AN][h_AN][i][j]; /* OptResidualH */
675 my_sum += r*orx;
676
677 }
678 }
679 }
680 }
681
682 } /* if (SpinP_switch==3 */
683
684 else{
685
686 for (spin=0; spin<=SpinP_switch; spin++){
687 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
688 Gc_AN = M2G[Mc_AN];
689 Cwan = WhatSpecies[Gc_AN];
690 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
691 Gh_AN = natn[Gc_AN][h_AN];
692 Hwan = WhatSpecies[Gh_AN];
693 for (i=0; i<Spe_Total_NO[Cwan]; i++){
694 for (j=0; j<Spe_Total_NO[Hwan]; j++){
695
696 s = HisH1[m][spin][Mc_AN][h_AN][i][j] - HisH1[m+1][spin][Mc_AN][h_AN][i][j];
697 y = ResidualH1[m][spin][Mc_AN][h_AN][i][j] - ResidualH1[m+1][spin][Mc_AN][h_AN][i][j];
698 r = s - al*y;
699 orx = ResidualH1[dim+1][spin][Mc_AN][h_AN][i][j];
700 my_sum += r*orx;
701 }
702 }
703 }
704 }
705 }
706
707 } /* else */
708
709 MPI_Allreduce(&my_sum, &ror[m], 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
710
711 } /* m */
712
713 /****************************************************
714 calculation of \sum_j b_{ij} * <r_j|OptResidualH>
715 ****************************************************/
716
717 for (m=0; m<dim; m++){
718 sum = 0.0;
719 for (n=0; n<dim; n++){
720 sum += IA[m][n]*ror[n];
721 }
722
723 coes2[m] = sum;
724 }
725
726 /****************************************************
727 mixing of Hamiltonian
728 ****************************************************/
729
730 if (1.0e-1<=NormRD[0])
731 alpha = 0.5;
732 else if (1.0e-2<=NormRD[0] && NormRD[0]<1.0e-1)
733 alpha = 0.6;
734 else if (1.0e-3<=NormRD[0] && NormRD[0]<1.0e-2)
735 alpha = 0.7;
736 else if (1.0e-4<=NormRD[0] && NormRD[0]<1.0e-3)
737 alpha = 0.8;
738 else
739 alpha = 1.0;
740
741 if (SpinP_switch==3){
742
743 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
744 Gc_AN = M2G[Mc_AN];
745 Cwan = WhatSpecies[Gc_AN];
746 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
747 Gh_AN = natn[Gc_AN][h_AN];
748 Hwan = WhatSpecies[Gh_AN];
749 for (i=0; i<Spe_Total_NO[Cwan]; i++){
750 for (j=0; j<Spe_Total_NO[Hwan]; j++){
751
752 h10 = 0.0;
753 h11 = 0.0;
754 h12 = 0.0;
755 h13 = 0.0;
756
757 h20 = 0.0;
758 h21 = 0.0;
759 h22 = 0.0;
760
761 for (m=0; m<dim; m++){
762
763 h10 += HisH1[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
764 h11 += HisH1[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
765 h12 += HisH1[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
766 h13 += HisH1[m][3][Mc_AN][h_AN][i][j]*coes[m+1];
767
768 h20 += HisH2[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
769 h21 += HisH2[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
770 h22 += HisH2[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
771
772 s = HisH1[m][0][Mc_AN][h_AN][i][j] - HisH1[m+1][0][Mc_AN][h_AN][i][j]; /* s */
773 y = ResidualH1[m][0][Mc_AN][h_AN][i][j] - ResidualH1[m+1][0][Mc_AN][h_AN][i][j]; /* y */
774 r = s - al*y; /* r */
775 h10 -= r*coes2[m];
776
777 s = HisH1[m][1][Mc_AN][h_AN][i][j] - HisH1[m+1][1][Mc_AN][h_AN][i][j]; /* s */
778 y = ResidualH1[m][1][Mc_AN][h_AN][i][j] - ResidualH1[m+1][1][Mc_AN][h_AN][i][j]; /* y */
779 r = s - al*y; /* r */
780 h11 -= r*coes2[m];
781
782 s = HisH1[m][2][Mc_AN][h_AN][i][j] - HisH1[m+1][2][Mc_AN][h_AN][i][j]; /* s */
783 y = ResidualH1[m][2][Mc_AN][h_AN][i][j] - ResidualH1[m+1][2][Mc_AN][h_AN][i][j]; /* y */
784 r = s - al*y; /* r */
785 h12 -= r*coes2[m];
786
787 s = HisH1[m][3][Mc_AN][h_AN][i][j] - HisH1[m+1][3][Mc_AN][h_AN][i][j]; /* s */
788 y = ResidualH1[m][3][Mc_AN][h_AN][i][j] - ResidualH1[m+1][3][Mc_AN][h_AN][i][j]; /* y */
789 r = s - al*y; /* r */
790 h13 -= r*coes2[m];
791
792 s = HisH2[m][0][Mc_AN][h_AN][i][j] - HisH2[m+1][0][Mc_AN][h_AN][i][j]; /* s */
793 y = ResidualH2[m][0][Mc_AN][h_AN][i][j] - ResidualH2[m+1][0][Mc_AN][h_AN][i][j]; /* y */
794 r = s - al*y; /* r */
795 h20 -= r*coes2[m];
796
797 s = HisH2[m][1][Mc_AN][h_AN][i][j] - HisH2[m+1][1][Mc_AN][h_AN][i][j]; /* s */
798 y = ResidualH2[m][1][Mc_AN][h_AN][i][j] - ResidualH2[m+1][1][Mc_AN][h_AN][i][j]; /* y */
799 r = s - al*y; /* r */
800 h21 -= r*coes2[m];
801
802 s = HisH2[m][2][Mc_AN][h_AN][i][j] - HisH2[m+1][2][Mc_AN][h_AN][i][j]; /* s */
803 y = ResidualH2[m][2][Mc_AN][h_AN][i][j] - ResidualH2[m+1][2][Mc_AN][h_AN][i][j]; /* y */
804 r = s - al*y; /* r */
805 h22 -= r*coes2[m];
806 }
807
808 H[0][Mc_AN][h_AN][i][j] = h10 - al*ResidualH1[dim+1][0][Mc_AN][h_AN][i][j];
809 H[1][Mc_AN][h_AN][i][j] = h11 - al*ResidualH1[dim+1][1][Mc_AN][h_AN][i][j];
810 H[2][Mc_AN][h_AN][i][j] = h12 - al*ResidualH1[dim+1][2][Mc_AN][h_AN][i][j];
811 H[3][Mc_AN][h_AN][i][j] = h13 - al*ResidualH1[dim+1][3][Mc_AN][h_AN][i][j];
812
813 iHNL[0][Mc_AN][h_AN][i][j] = h20 - al*ResidualH2[dim+1][0][Mc_AN][h_AN][i][j];
814 iHNL[1][Mc_AN][h_AN][i][j] = h21 - al*ResidualH2[dim+1][1][Mc_AN][h_AN][i][j];
815 iHNL[2][Mc_AN][h_AN][i][j] = h22 - al*ResidualH2[dim+1][2][Mc_AN][h_AN][i][j];
816
817 }
818 }
819 }
820 }
821
822 }
823
824 else{
825
826 for (spin=0; spin<=SpinP_switch; spin++){
827 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
828 Gc_AN = M2G[Mc_AN];
829 Cwan = WhatSpecies[Gc_AN];
830 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
831 Gh_AN = natn[Gc_AN][h_AN];
832 Hwan = WhatSpecies[Gh_AN];
833 for (i=0; i<Spe_Total_NO[Cwan]; i++){
834 for (j=0; j<Spe_Total_NO[Hwan]; j++){
835
836 h = 0.0;
837 for (m=0; m<dim; m++){
838
839 h += HisH1[m][spin][Mc_AN][h_AN][i][j]*coes[m+1];
840 s = HisH1[m][spin][Mc_AN][h_AN][i][j] - HisH1[m+1][spin][Mc_AN][h_AN][i][j];
841 y = ResidualH1[m][spin][Mc_AN][h_AN][i][j] - ResidualH1[m+1][spin][Mc_AN][h_AN][i][j];
842 r = s - al*y;
843 h -= r*coes2[m];
844 }
845
846 H[spin][Mc_AN][h_AN][i][j] = h - al*ResidualH1[dim+1][spin][Mc_AN][h_AN][i][j];
847 }
848 }
849 }
850 }
851 }
852 }
853
854 /****************************************************
855 shifting of Hamiltonian
856 ****************************************************/
857
858 if (SpinP_switch==3){
859
860 /* shift the current Hamiltonian */
861
862 for (m=dim; 0<m; m--){
863 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
864 Gc_AN = M2G[Mc_AN];
865 Cwan = WhatSpecies[Gc_AN];
866 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
867 Gh_AN = natn[Gc_AN][h_AN];
868 Hwan = WhatSpecies[Gh_AN];
869 for (i=0; i<Spe_Total_NO[Cwan]; i++){
870 for (j=0; j<Spe_Total_NO[Hwan]; j++){
871
872 HisH1[m][0][Mc_AN][h_AN][i][j] = HisH1[m-1][0][Mc_AN][h_AN][i][j];
873 HisH1[m][1][Mc_AN][h_AN][i][j] = HisH1[m-1][1][Mc_AN][h_AN][i][j];
874 HisH1[m][2][Mc_AN][h_AN][i][j] = HisH1[m-1][2][Mc_AN][h_AN][i][j];
875 HisH1[m][3][Mc_AN][h_AN][i][j] = HisH1[m-1][3][Mc_AN][h_AN][i][j];
876
877 HisH2[m][0][Mc_AN][h_AN][i][j] = HisH2[m-1][0][Mc_AN][h_AN][i][j];
878 HisH2[m][1][Mc_AN][h_AN][i][j] = HisH2[m-1][1][Mc_AN][h_AN][i][j];
879 HisH2[m][2][Mc_AN][h_AN][i][j] = HisH2[m-1][2][Mc_AN][h_AN][i][j];
880 }
881 }
882 }
883 }
884 }
885
886 /* save the current Hamiltonian */
887
888 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
889 Gc_AN = M2G[Mc_AN];
890 Cwan = WhatSpecies[Gc_AN];
891 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
892 Gh_AN = natn[Gc_AN][h_AN];
893 Hwan = WhatSpecies[Gh_AN];
894 for (i=0; i<Spe_Total_NO[Cwan]; i++){
895 for (j=0; j<Spe_Total_NO[Hwan]; j++){
896
897 HisH1[0][0][Mc_AN][h_AN][i][j] = H[0][Mc_AN][h_AN][i][j];
898 HisH1[0][1][Mc_AN][h_AN][i][j] = H[1][Mc_AN][h_AN][i][j];
899 HisH1[0][2][Mc_AN][h_AN][i][j] = H[2][Mc_AN][h_AN][i][j];
900 HisH1[0][3][Mc_AN][h_AN][i][j] = H[3][Mc_AN][h_AN][i][j];
901
902 HisH2[0][0][Mc_AN][h_AN][i][j] = iHNL[0][Mc_AN][h_AN][i][j];
903 HisH2[0][1][Mc_AN][h_AN][i][j] = iHNL[1][Mc_AN][h_AN][i][j];
904 HisH2[0][2][Mc_AN][h_AN][i][j] = iHNL[2][Mc_AN][h_AN][i][j];
905 }
906 }
907 }
908 }
909
910 }
911
912 else {
913
914 /* shift the current Hamiltonian */
915
916 for (m=dim; 0<m; m--){
917 for (spin=0; spin<=SpinP_switch; spin++){
918 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
919 Gc_AN = M2G[Mc_AN];
920 Cwan = WhatSpecies[Gc_AN];
921 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
922 Gh_AN = natn[Gc_AN][h_AN];
923 Hwan = WhatSpecies[Gh_AN];
924 for (i=0; i<Spe_Total_NO[Cwan]; i++){
925 for (j=0; j<Spe_Total_NO[Hwan]; j++){
926 HisH1[m][spin][Mc_AN][h_AN][i][j] = HisH1[m-1][spin][Mc_AN][h_AN][i][j];
927 }
928 }
929 }
930 }
931 }
932 }
933
934 /* save the current Hamiltonian */
935
936 for (spin=0; spin<=SpinP_switch; spin++){
937 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
938 Gc_AN = M2G[Mc_AN];
939 Cwan = WhatSpecies[Gc_AN];
940 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
941 Gh_AN = natn[Gc_AN][h_AN];
942 Hwan = WhatSpecies[Gh_AN];
943 for (i=0; i<Spe_Total_NO[Cwan]; i++){
944 for (j=0; j<Spe_Total_NO[Hwan]; j++){
945 HisH1[0][spin][Mc_AN][h_AN][i][j] = H[spin][Mc_AN][h_AN][i][j];
946 }
947 }
948 }
949 }
950 }
951
952 } /* else */
953
954 /****************************************************
955 freeing of arrays
956 ****************************************************/
957
958 free(coes);
959 free(coes2);
960 free(ror);
961
962 for (i=0; i<List_YOUSO[39]; i++){
963 free(A[i]);
964 }
965 free(A);
966
967 for (i=0; i<List_YOUSO[39]; i++){
968 free(IA[i]);
969 }
970 free(IA);
971 }
972
973
974
975
976
977
978
979
980
981
Pulay_Mixing_H_with_One_Shot_Hessian(int MD_iter,int SCF_iter,int SCF_iter0)982 void Pulay_Mixing_H_with_One_Shot_Hessian(int MD_iter, int SCF_iter, int SCF_iter0 )
983 {
984 int Mc_AN,Gc_AN,Cwan,Hwan,h_AN,Gh_AN,i,j,spin;
985 int dim,m,n,flag_nan;
986 double my_sum,tmp1,tmp2,alpha;
987 double r,r10,r11,r12,r13,r20,r21,r22;
988 double h,h10,h11,h12,h13,h20,h21,h22;
989 double my_sy,my_yy,sy,yy,norm,s,y,orx,al,be;
990 double **A,**IA,*coes;
991 char nanchar[300];
992
993 /****************************************************
994 determination of dimension of the subspace
995 ****************************************************/
996
997 if (SCF_iter<=Num_Mixing_pDM) dim = SCF_iter-1;
998 else dim = Num_Mixing_pDM;
999
1000 /****************************************************
1001 allocation of arrays
1002 ****************************************************/
1003
1004 coes = (double*)malloc(sizeof(double)*List_YOUSO[39]);
1005
1006 A = (double**)malloc(sizeof(double*)*List_YOUSO[39]);
1007 for (i=0; i<List_YOUSO[39]; i++){
1008 A[i] = (double*)malloc(sizeof(double)*List_YOUSO[39]);
1009 }
1010
1011 IA = (double**)malloc(sizeof(double*)*List_YOUSO[39]);
1012 for (i=0; i<List_YOUSO[39]; i++){
1013 IA[i] = (double*)malloc(sizeof(double)*List_YOUSO[39]);
1014 }
1015
1016 /****************************************************
1017 shift the residual H
1018 ****************************************************/
1019
1020 if (SpinP_switch==3){
1021
1022 /* shift the residual Hamiltonian */
1023
1024 for (m=dim; 0<m; m--){
1025 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1026 Gc_AN = M2G[Mc_AN];
1027 Cwan = WhatSpecies[Gc_AN];
1028 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1029 Gh_AN = natn[Gc_AN][h_AN];
1030 Hwan = WhatSpecies[Gh_AN];
1031 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1032 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1033
1034 ResidualH1[m][0][Mc_AN][h_AN][i][j] = ResidualH1[m-1][0][Mc_AN][h_AN][i][j];
1035 ResidualH1[m][1][Mc_AN][h_AN][i][j] = ResidualH1[m-1][1][Mc_AN][h_AN][i][j];
1036 ResidualH1[m][2][Mc_AN][h_AN][i][j] = ResidualH1[m-1][2][Mc_AN][h_AN][i][j];
1037 ResidualH1[m][3][Mc_AN][h_AN][i][j] = ResidualH1[m-1][3][Mc_AN][h_AN][i][j];
1038
1039 ResidualH2[m][0][Mc_AN][h_AN][i][j] = ResidualH2[m-1][0][Mc_AN][h_AN][i][j];
1040 ResidualH2[m][1][Mc_AN][h_AN][i][j] = ResidualH2[m-1][1][Mc_AN][h_AN][i][j];
1041 ResidualH2[m][2][Mc_AN][h_AN][i][j] = ResidualH2[m-1][2][Mc_AN][h_AN][i][j];
1042 }
1043 }
1044 }
1045 }
1046 }
1047
1048 /* calculate the current residual Hamiltonian */
1049
1050 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1051 Gc_AN = M2G[Mc_AN];
1052 Cwan = WhatSpecies[Gc_AN];
1053 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1054 Gh_AN = natn[Gc_AN][h_AN];
1055 Hwan = WhatSpecies[Gh_AN];
1056 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1057 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1058
1059 ResidualH1[0][0][Mc_AN][h_AN][i][j] = H[0][Mc_AN][h_AN][i][j] - HisH1[0][0][Mc_AN][h_AN][i][j];
1060 ResidualH1[0][1][Mc_AN][h_AN][i][j] = H[1][Mc_AN][h_AN][i][j] - HisH1[0][1][Mc_AN][h_AN][i][j];
1061 ResidualH1[0][2][Mc_AN][h_AN][i][j] = H[2][Mc_AN][h_AN][i][j] - HisH1[0][2][Mc_AN][h_AN][i][j];
1062 ResidualH1[0][3][Mc_AN][h_AN][i][j] = H[3][Mc_AN][h_AN][i][j] - HisH1[0][3][Mc_AN][h_AN][i][j];
1063
1064 ResidualH2[0][0][Mc_AN][h_AN][i][j] = iHNL[0][Mc_AN][h_AN][i][j] - HisH2[0][0][Mc_AN][h_AN][i][j];
1065 ResidualH2[0][1][Mc_AN][h_AN][i][j] = iHNL[1][Mc_AN][h_AN][i][j] - HisH2[0][1][Mc_AN][h_AN][i][j];
1066 ResidualH2[0][2][Mc_AN][h_AN][i][j] = iHNL[2][Mc_AN][h_AN][i][j] - HisH2[0][2][Mc_AN][h_AN][i][j];
1067 }
1068 }
1069 }
1070 }
1071
1072 }
1073
1074 else{
1075
1076 /* shift the residual Hamiltonian */
1077
1078 for (m=dim; 0<m; m--){
1079 for (spin=0; spin<=SpinP_switch; spin++){
1080 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1081 Gc_AN = M2G[Mc_AN];
1082 Cwan = WhatSpecies[Gc_AN];
1083 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1084 Gh_AN = natn[Gc_AN][h_AN];
1085 Hwan = WhatSpecies[Gh_AN];
1086 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1087 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1088 ResidualH1[m][spin][Mc_AN][h_AN][i][j] = ResidualH1[m-1][spin][Mc_AN][h_AN][i][j];
1089 }
1090 }
1091 }
1092 }
1093 }
1094 }
1095
1096 /* calculate the current residual Hamiltonian */
1097
1098 for (spin=0; spin<=SpinP_switch; spin++){
1099 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1100 Gc_AN = M2G[Mc_AN];
1101 Cwan = WhatSpecies[Gc_AN];
1102 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1103 Gh_AN = natn[Gc_AN][h_AN];
1104 Hwan = WhatSpecies[Gh_AN];
1105 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1106 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1107 ResidualH1[0][spin][Mc_AN][h_AN][i][j] = H[spin][Mc_AN][h_AN][i][j] - HisH1[0][spin][Mc_AN][h_AN][i][j];
1108 }
1109 }
1110 }
1111 }
1112 }
1113
1114 } /* else */
1115
1116 /****************************************************
1117 calculation of the residual matrix
1118 ****************************************************/
1119
1120 for (m=0; m<dim; m++){
1121 for (n=0; n<dim; n++){
1122
1123 my_sum = 0.0;
1124
1125 if (SpinP_switch==3){
1126
1127 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1128 Gc_AN = M2G[Mc_AN];
1129 Cwan = WhatSpecies[Gc_AN];
1130 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1131 Gh_AN = natn[Gc_AN][h_AN];
1132 Hwan = WhatSpecies[Gh_AN];
1133 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1134 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1135
1136 tmp1 = ResidualH1[m][0][Mc_AN][h_AN][i][j];
1137 tmp2 = ResidualH1[n][0][Mc_AN][h_AN][i][j];
1138 my_sum += tmp1*tmp2;
1139
1140 tmp1 = ResidualH1[m][1][Mc_AN][h_AN][i][j];
1141 tmp2 = ResidualH1[n][1][Mc_AN][h_AN][i][j];
1142 my_sum += tmp1*tmp2;
1143
1144 tmp1 = ResidualH1[m][2][Mc_AN][h_AN][i][j];
1145 tmp2 = ResidualH1[n][2][Mc_AN][h_AN][i][j];
1146 my_sum += tmp1*tmp2;
1147
1148 tmp1 = ResidualH1[m][3][Mc_AN][h_AN][i][j];
1149 tmp2 = ResidualH1[n][3][Mc_AN][h_AN][i][j];
1150 my_sum += tmp1*tmp2;
1151
1152 tmp1 = ResidualH2[m][0][Mc_AN][h_AN][i][j];
1153 tmp2 = ResidualH2[n][0][Mc_AN][h_AN][i][j];
1154 my_sum += tmp1*tmp2;
1155
1156 tmp1 = ResidualH2[m][1][Mc_AN][h_AN][i][j];
1157 tmp2 = ResidualH2[n][1][Mc_AN][h_AN][i][j];
1158 my_sum += tmp1*tmp2;
1159
1160 tmp1 = ResidualH2[m][2][Mc_AN][h_AN][i][j];
1161 tmp2 = ResidualH2[n][2][Mc_AN][h_AN][i][j];
1162 my_sum += tmp1*tmp2;
1163 }
1164 }
1165 }
1166 }
1167
1168 } /* if (SpinP_switch==3 */
1169
1170 else{
1171
1172 for (spin=0; spin<=SpinP_switch; spin++){
1173 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1174 Gc_AN = M2G[Mc_AN];
1175 Cwan = WhatSpecies[Gc_AN];
1176 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1177 Gh_AN = natn[Gc_AN][h_AN];
1178 Hwan = WhatSpecies[Gh_AN];
1179 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1180 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1181 tmp1 = ResidualH1[m][spin][Mc_AN][h_AN][i][j];
1182 tmp2 = ResidualH1[n][spin][Mc_AN][h_AN][i][j];
1183 my_sum += tmp1*tmp2;
1184 }
1185 }
1186 }
1187 }
1188 }
1189
1190 } /* else */
1191
1192 MPI_Allreduce(&my_sum, &A[m][n], 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
1193 A[n][m] = A[m][n];
1194
1195 } /* n */
1196 } /* m */
1197
1198 NormRD[0] = A[0][0];
1199
1200 for (m=1; m<=dim; m++){
1201 A[m-1][dim] = -1.0;
1202 A[dim][m-1] = -1.0;
1203 }
1204 A[dim][dim] = 0.0;
1205
1206 Inverse(dim,A,IA);
1207
1208 for (m=1; m<=dim; m++){
1209 coes[m] = -IA[m-1][dim];
1210 }
1211
1212 /****************************************************
1213 check "nan", "NaN", "inf" or "Inf"
1214 ****************************************************/
1215
1216 flag_nan = 0;
1217 for (m=1; m<=dim; m++){
1218
1219 sprintf(nanchar,"%8.4f",coes[m]);
1220 if ( strstr(nanchar,"nan")!=NULL || strstr(nanchar,"NaN")!=NULL
1221 || strstr(nanchar,"inf")!=NULL || strstr(nanchar,"Inf")!=NULL){
1222
1223 flag_nan = 1;
1224 }
1225 }
1226
1227 if (flag_nan==1){
1228 for (m=1; m<=dim; m++){
1229 coes[m] = 0.0;
1230 }
1231 coes[1] = 0.05;
1232 coes[2] = 0.95;
1233 }
1234
1235 /****************************************************
1236 calculation of optimum residual Hamiltonian
1237 ****************************************************/
1238
1239 if (SpinP_switch==3){
1240
1241 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1242 Gc_AN = M2G[Mc_AN];
1243 Cwan = WhatSpecies[Gc_AN];
1244 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1245 Gh_AN = natn[Gc_AN][h_AN];
1246 Hwan = WhatSpecies[Gh_AN];
1247 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1248 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1249
1250 r10 = 0.0;
1251 r11 = 0.0;
1252 r12 = 0.0;
1253 r13 = 0.0;
1254
1255 r20 = 0.0;
1256 r21 = 0.0;
1257 r22 = 0.0;
1258
1259 for (m=0; m<dim; m++){
1260
1261 r10 += ResidualH1[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
1262 r11 += ResidualH1[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
1263 r12 += ResidualH1[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
1264 r13 += ResidualH1[m][3][Mc_AN][h_AN][i][j]*coes[m+1];
1265
1266 r20 += ResidualH2[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
1267 r21 += ResidualH2[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
1268 r22 += ResidualH2[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
1269 }
1270
1271 /* optimum Residual H is stored in ResidualH1[dim] and ResidualH2[dim] */
1272
1273 ResidualH1[dim][0][Mc_AN][h_AN][i][j] = r10;
1274 ResidualH1[dim][1][Mc_AN][h_AN][i][j] = r11;
1275 ResidualH1[dim][2][Mc_AN][h_AN][i][j] = r12;
1276 ResidualH1[dim][3][Mc_AN][h_AN][i][j] = r13;
1277
1278 ResidualH2[dim][0][Mc_AN][h_AN][i][j] = r20;
1279 ResidualH2[dim][1][Mc_AN][h_AN][i][j] = r21;
1280 ResidualH2[dim][2][Mc_AN][h_AN][i][j] = r22;
1281
1282 }
1283 }
1284 }
1285 }
1286
1287 }
1288
1289 else{
1290
1291 for (spin=0; spin<=SpinP_switch; spin++){
1292 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1293 Gc_AN = M2G[Mc_AN];
1294 Cwan = WhatSpecies[Gc_AN];
1295 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1296 Gh_AN = natn[Gc_AN][h_AN];
1297 Hwan = WhatSpecies[Gh_AN];
1298 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1299 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1300
1301 r = 0.0;
1302 for (m=0; m<dim; m++){
1303 r += ResidualH1[m][spin][Mc_AN][h_AN][i][j]*coes[m+1];
1304 }
1305
1306 /* optimum Residual H is stored in ResidualH1[dim] */
1307
1308 ResidualH1[dim][spin][Mc_AN][h_AN][i][j] = r;
1309
1310 }
1311 }
1312 }
1313 }
1314 }
1315 }
1316
1317 /****************************************************
1318 innner products of <s|y> and <y|y>
1319 ****************************************************/
1320
1321 my_sy = 0.0;
1322 my_yy = 0.0;
1323
1324 if (SpinP_switch==3){
1325
1326 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1327 Gc_AN = M2G[Mc_AN];
1328 Cwan = WhatSpecies[Gc_AN];
1329 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1330 Gh_AN = natn[Gc_AN][h_AN];
1331 Hwan = WhatSpecies[Gh_AN];
1332 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1333 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1334
1335 tmp1 = HisH1[0][0][Mc_AN][h_AN][i][j] - HisH1[1][0][Mc_AN][h_AN][i][j]; /* s */
1336 tmp2 = ResidualH1[0][0][Mc_AN][h_AN][i][j] - ResidualH1[1][0][Mc_AN][h_AN][i][j]; /* y */
1337 my_sy += tmp1*tmp2;
1338 my_yy += tmp2*tmp2;
1339
1340 tmp1 = HisH1[0][1][Mc_AN][h_AN][i][j] - HisH1[1][1][Mc_AN][h_AN][i][j]; /* s */
1341 tmp2 = ResidualH1[0][1][Mc_AN][h_AN][i][j] - ResidualH1[1][1][Mc_AN][h_AN][i][j]; /* y */
1342 my_sy += tmp1*tmp2;
1343 my_yy += tmp2*tmp2;
1344
1345 tmp1 = HisH1[0][2][Mc_AN][h_AN][i][j] - HisH1[1][2][Mc_AN][h_AN][i][j]; /* s */
1346 tmp2 = ResidualH1[0][2][Mc_AN][h_AN][i][j] - ResidualH1[1][2][Mc_AN][h_AN][i][j]; /* y */
1347 my_sy += tmp1*tmp2;
1348 my_yy += tmp2*tmp2;
1349
1350 tmp1 = HisH1[0][3][Mc_AN][h_AN][i][j] - HisH1[1][3][Mc_AN][h_AN][i][j]; /* s */
1351 tmp2 = ResidualH1[0][3][Mc_AN][h_AN][i][j] - ResidualH1[1][3][Mc_AN][h_AN][i][j]; /* y */
1352 my_sy += tmp1*tmp2;
1353 my_yy += tmp2*tmp2;
1354
1355 tmp1 = HisH2[0][0][Mc_AN][h_AN][i][j] - HisH2[1][0][Mc_AN][h_AN][i][j]; /* s */
1356 tmp2 = ResidualH2[0][0][Mc_AN][h_AN][i][j] - ResidualH2[1][0][Mc_AN][h_AN][i][j]; /* y */
1357 my_sy += tmp1*tmp2;
1358 my_yy += tmp2*tmp2;
1359
1360 tmp1 = HisH2[0][1][Mc_AN][h_AN][i][j] - HisH2[1][1][Mc_AN][h_AN][i][j]; /* s */
1361 tmp2 = ResidualH2[0][1][Mc_AN][h_AN][i][j] - ResidualH2[1][1][Mc_AN][h_AN][i][j]; /* y */
1362 my_sy += tmp1*tmp2;
1363 my_yy += tmp2*tmp2;
1364
1365 tmp1 = HisH2[0][2][Mc_AN][h_AN][i][j] - HisH2[1][2][Mc_AN][h_AN][i][j]; /* s */
1366 tmp2 = ResidualH2[0][2][Mc_AN][h_AN][i][j] - ResidualH2[1][2][Mc_AN][h_AN][i][j]; /* y */
1367 my_sy += tmp1*tmp2;
1368 my_yy += tmp2*tmp2;
1369 }
1370 }
1371 }
1372 }
1373
1374 } /* if (SpinP_switch==3 */
1375
1376 else{
1377
1378 for (spin=0; spin<=SpinP_switch; spin++){
1379 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1380 Gc_AN = M2G[Mc_AN];
1381 Cwan = WhatSpecies[Gc_AN];
1382 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1383 Gh_AN = natn[Gc_AN][h_AN];
1384 Hwan = WhatSpecies[Gh_AN];
1385 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1386 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1387
1388 tmp1 = HisH1[0][spin][Mc_AN][h_AN][i][j] - HisH1[1][spin][Mc_AN][h_AN][i][j]; /* s */
1389 tmp2 = ResidualH1[0][spin][Mc_AN][h_AN][i][j] - ResidualH1[1][spin][Mc_AN][h_AN][i][j]; /* y */
1390 my_sy += tmp1*tmp2;
1391 my_yy += tmp2*tmp2;
1392 }
1393 }
1394 }
1395 }
1396 }
1397
1398 } /* else */
1399
1400 MPI_Allreduce(&my_sy, &sy, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
1401 MPI_Allreduce(&my_yy, &yy, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
1402
1403 /* al < sy/yy */
1404
1405 al = sy/yy - 0.2;
1406
1407 /* be = 1/(<s|y>-al*<y|y>) */
1408
1409 be = 1.0/(sy-al*yy);
1410
1411 /****************************************************
1412 inner product of (<s|-al<y|)|OptResidualH>
1413 ****************************************************/
1414
1415 my_sum = 0.0;
1416
1417 if (SpinP_switch==3){
1418
1419 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1420 Gc_AN = M2G[Mc_AN];
1421 Cwan = WhatSpecies[Gc_AN];
1422 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1423 Gh_AN = natn[Gc_AN][h_AN];
1424 Hwan = WhatSpecies[Gh_AN];
1425 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1426 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1427
1428 s = HisH1[0][0][Mc_AN][h_AN][i][j] - HisH1[1][0][Mc_AN][h_AN][i][j]; /* s */
1429 y = ResidualH1[0][0][Mc_AN][h_AN][i][j] - ResidualH1[1][0][Mc_AN][h_AN][i][j]; /* y */
1430 orx = ResidualH1[dim][0][Mc_AN][h_AN][i][j]; /* OptResidualH */
1431 my_sum += (s-al*y)*orx;
1432
1433 s = HisH1[0][1][Mc_AN][h_AN][i][j] - HisH1[1][1][Mc_AN][h_AN][i][j]; /* s */
1434 y = ResidualH1[0][1][Mc_AN][h_AN][i][j] - ResidualH1[1][1][Mc_AN][h_AN][i][j]; /* y */
1435 orx = ResidualH1[dim][1][Mc_AN][h_AN][i][j]; /* OptResidualH */
1436 my_sum += (s-al*y)*orx;
1437
1438 s = HisH1[0][2][Mc_AN][h_AN][i][j] - HisH1[1][2][Mc_AN][h_AN][i][j]; /* s */
1439 y = ResidualH1[0][2][Mc_AN][h_AN][i][j] - ResidualH1[1][2][Mc_AN][h_AN][i][j]; /* y */
1440 orx = ResidualH1[dim][2][Mc_AN][h_AN][i][j]; /* OptResidualH */
1441 my_sum += (s-al*y)*orx;
1442
1443 s = HisH1[0][3][Mc_AN][h_AN][i][j] - HisH1[1][3][Mc_AN][h_AN][i][j]; /* s */
1444 y = ResidualH1[0][3][Mc_AN][h_AN][i][j] - ResidualH1[1][3][Mc_AN][h_AN][i][j]; /* y */
1445 orx = ResidualH1[dim][3][Mc_AN][h_AN][i][j]; /* OptResidualH */
1446 my_sum += (s-al*y)*orx;
1447
1448 s = HisH2[0][0][Mc_AN][h_AN][i][j] - HisH2[1][0][Mc_AN][h_AN][i][j]; /* s */
1449 y = ResidualH2[0][0][Mc_AN][h_AN][i][j] - ResidualH2[1][0][Mc_AN][h_AN][i][j]; /* y */
1450 orx = ResidualH2[dim][0][Mc_AN][h_AN][i][j]; /* OptResidualH */
1451 my_sum += (s-al*y)*orx;
1452
1453 s = HisH2[0][1][Mc_AN][h_AN][i][j] - HisH2[1][1][Mc_AN][h_AN][i][j]; /* s */
1454 y = ResidualH2[0][1][Mc_AN][h_AN][i][j] - ResidualH2[1][1][Mc_AN][h_AN][i][j]; /* y */
1455 orx = ResidualH2[dim][1][Mc_AN][h_AN][i][j]; /* OptResidualH */
1456 my_sum += (s-al*y)*orx;
1457
1458 s = HisH2[0][2][Mc_AN][h_AN][i][j] - HisH2[1][2][Mc_AN][h_AN][i][j]; /* s */
1459 y = ResidualH2[0][2][Mc_AN][h_AN][i][j] - ResidualH2[1][2][Mc_AN][h_AN][i][j]; /* y */
1460 orx = ResidualH2[dim][2][Mc_AN][h_AN][i][j]; /* OptResidualH */
1461 my_sum += (s-al*y)*orx;
1462 }
1463 }
1464 }
1465 }
1466
1467 } /* if (SpinP_switch==3 */
1468
1469 else{
1470
1471 for (spin=0; spin<=SpinP_switch; spin++){
1472 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1473 Gc_AN = M2G[Mc_AN];
1474 Cwan = WhatSpecies[Gc_AN];
1475 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1476 Gh_AN = natn[Gc_AN][h_AN];
1477 Hwan = WhatSpecies[Gh_AN];
1478 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1479 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1480 s = HisH1[0][spin][Mc_AN][h_AN][i][j] - HisH1[1][spin][Mc_AN][h_AN][i][j]; /* s */
1481 y = ResidualH1[0][spin][Mc_AN][h_AN][i][j] - ResidualH1[1][spin][Mc_AN][h_AN][i][j]; /* y */
1482 orx = ResidualH1[dim][spin][Mc_AN][h_AN][i][j]; /* OptResidualH */
1483 my_sum += (s-al*y)*orx;
1484 }
1485 }
1486 }
1487 }
1488 }
1489
1490 } /* else */
1491
1492 MPI_Allreduce(&my_sum, &norm, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
1493 be = norm*be;
1494
1495 /****************************************************
1496 mixing of Hamiltonian
1497 ****************************************************/
1498
1499 if (1.0e-1<=NormRD[0])
1500 alpha = 0.5;
1501 else if (1.0e-2<=NormRD[0] && NormRD[0]<1.0e-1)
1502 alpha = 0.6;
1503 else if (1.0e-3<=NormRD[0] && NormRD[0]<1.0e-2)
1504 alpha = 0.7;
1505 else if (1.0e-4<=NormRD[0] && NormRD[0]<1.0e-3)
1506 alpha = 0.8;
1507 else
1508 alpha = 1.0;
1509
1510 if (SpinP_switch==3){
1511
1512 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1513 Gc_AN = M2G[Mc_AN];
1514 Cwan = WhatSpecies[Gc_AN];
1515 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1516 Gh_AN = natn[Gc_AN][h_AN];
1517 Hwan = WhatSpecies[Gh_AN];
1518 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1519 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1520
1521 h10 = 0.0;
1522 h11 = 0.0;
1523 h12 = 0.0;
1524 h13 = 0.0;
1525
1526 h20 = 0.0;
1527 h21 = 0.0;
1528 h22 = 0.0;
1529
1530 for (m=0; m<dim; m++){
1531
1532 h10 += HisH1[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
1533 h11 += HisH1[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
1534 h12 += HisH1[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
1535 h13 += HisH1[m][3][Mc_AN][h_AN][i][j]*coes[m+1];
1536
1537 h20 += HisH2[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
1538 h21 += HisH2[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
1539 h22 += HisH2[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
1540 }
1541
1542 s = HisH1[0][0][Mc_AN][h_AN][i][j] - HisH1[1][0][Mc_AN][h_AN][i][j]; /* s */
1543 y = ResidualH1[0][0][Mc_AN][h_AN][i][j] - ResidualH1[1][0][Mc_AN][h_AN][i][j]; /* y */
1544 orx = ResidualH1[dim][0][Mc_AN][h_AN][i][j]; /* OptResidualH */
1545 H[0][Mc_AN][h_AN][i][j] = h10 - alpha*(al*orx + (s-al*y)*be);
1546
1547 s = HisH1[0][1][Mc_AN][h_AN][i][j] - HisH1[1][1][Mc_AN][h_AN][i][j]; /* s */
1548 y = ResidualH1[0][1][Mc_AN][h_AN][i][j] - ResidualH1[1][1][Mc_AN][h_AN][i][j]; /* y */
1549 orx = ResidualH1[dim][1][Mc_AN][h_AN][i][j]; /* OptResidualH */
1550 H[1][Mc_AN][h_AN][i][j] = h11 - alpha*(al*orx + (s-al*y)*be);
1551
1552 s = HisH1[0][2][Mc_AN][h_AN][i][j] - HisH1[1][2][Mc_AN][h_AN][i][j]; /* s */
1553 y = ResidualH1[0][2][Mc_AN][h_AN][i][j] - ResidualH1[1][2][Mc_AN][h_AN][i][j]; /* y */
1554 orx = ResidualH1[dim][2][Mc_AN][h_AN][i][j]; /* OptResidualH */
1555 H[2][Mc_AN][h_AN][i][j] = h12 - alpha*(al*orx + (s-al*y)*be);
1556
1557 s = HisH1[0][3][Mc_AN][h_AN][i][j] - HisH1[1][3][Mc_AN][h_AN][i][j]; /* s */
1558 y = ResidualH1[0][3][Mc_AN][h_AN][i][j] - ResidualH1[1][3][Mc_AN][h_AN][i][j]; /* y */
1559 orx = ResidualH1[dim][3][Mc_AN][h_AN][i][j]; /* OptResidualH */
1560 H[3][Mc_AN][h_AN][i][j] = h13 - alpha*(al*orx + (s-al*y)*be);
1561
1562 s = HisH2[0][0][Mc_AN][h_AN][i][j] - HisH2[1][0][Mc_AN][h_AN][i][j]; /* s */
1563 y = ResidualH2[0][0][Mc_AN][h_AN][i][j] - ResidualH2[1][0][Mc_AN][h_AN][i][j]; /* y */
1564 orx = ResidualH2[dim][0][Mc_AN][h_AN][i][j]; /* OptResidualH */
1565 iHNL[0][Mc_AN][h_AN][i][j] = h20 - alpha*(al*orx + (s-al*y)*be);
1566
1567 s = HisH2[0][1][Mc_AN][h_AN][i][j] - HisH2[1][1][Mc_AN][h_AN][i][j]; /* s */
1568 y = ResidualH2[0][1][Mc_AN][h_AN][i][j] - ResidualH2[1][1][Mc_AN][h_AN][i][j]; /* y */
1569 orx = ResidualH2[dim][1][Mc_AN][h_AN][i][j]; /* OptResidualH */
1570 iHNL[1][Mc_AN][h_AN][i][j] = h21 - alpha*(al*orx + (s-al*y)*be);
1571
1572 s = HisH2[0][2][Mc_AN][h_AN][i][j] - HisH2[1][2][Mc_AN][h_AN][i][j]; /* s */
1573 y = ResidualH2[0][2][Mc_AN][h_AN][i][j] - ResidualH2[1][2][Mc_AN][h_AN][i][j]; /* y */
1574 orx = ResidualH2[dim][2][Mc_AN][h_AN][i][j]; /* OptResidualH */
1575 iHNL[2][Mc_AN][h_AN][i][j] = h22 - alpha*(al*orx + (s-al*y)*be);
1576 }
1577 }
1578 }
1579 }
1580
1581 }
1582
1583 else{
1584
1585 for (spin=0; spin<=SpinP_switch; spin++){
1586 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1587 Gc_AN = M2G[Mc_AN];
1588 Cwan = WhatSpecies[Gc_AN];
1589 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1590 Gh_AN = natn[Gc_AN][h_AN];
1591 Hwan = WhatSpecies[Gh_AN];
1592 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1593 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1594
1595 h = 0.0;
1596 for (m=0; m<dim; m++){
1597 h += HisH1[m][spin][Mc_AN][h_AN][i][j]*coes[m+1];
1598 }
1599
1600 s = HisH1[0][spin][Mc_AN][h_AN][i][j] - HisH1[1][spin][Mc_AN][h_AN][i][j]; /* s */
1601 y = ResidualH1[0][spin][Mc_AN][h_AN][i][j] - ResidualH1[1][spin][Mc_AN][h_AN][i][j]; /* y */
1602 orx = ResidualH1[dim][spin][Mc_AN][h_AN][i][j]; /* OptResidualH */
1603 H[spin][Mc_AN][h_AN][i][j] = h - alpha*(al*orx + (s-al*y)*be);
1604
1605 }
1606 }
1607 }
1608 }
1609 }
1610 }
1611
1612 /****************************************************
1613 shifting of Hamiltonian
1614 ****************************************************/
1615
1616 if (SpinP_switch==3){
1617
1618 /* shift the current Hamiltonian */
1619
1620 for (m=dim; 0<m; m--){
1621 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1622 Gc_AN = M2G[Mc_AN];
1623 Cwan = WhatSpecies[Gc_AN];
1624 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1625 Gh_AN = natn[Gc_AN][h_AN];
1626 Hwan = WhatSpecies[Gh_AN];
1627 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1628 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1629
1630 HisH1[m][0][Mc_AN][h_AN][i][j] = HisH1[m-1][0][Mc_AN][h_AN][i][j];
1631 HisH1[m][1][Mc_AN][h_AN][i][j] = HisH1[m-1][1][Mc_AN][h_AN][i][j];
1632 HisH1[m][2][Mc_AN][h_AN][i][j] = HisH1[m-1][2][Mc_AN][h_AN][i][j];
1633 HisH1[m][3][Mc_AN][h_AN][i][j] = HisH1[m-1][3][Mc_AN][h_AN][i][j];
1634
1635 HisH2[m][0][Mc_AN][h_AN][i][j] = HisH2[m-1][0][Mc_AN][h_AN][i][j];
1636 HisH2[m][1][Mc_AN][h_AN][i][j] = HisH2[m-1][1][Mc_AN][h_AN][i][j];
1637 HisH2[m][2][Mc_AN][h_AN][i][j] = HisH2[m-1][2][Mc_AN][h_AN][i][j];
1638 }
1639 }
1640 }
1641 }
1642 }
1643
1644 /* save the current Hamiltonian */
1645
1646 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1647 Gc_AN = M2G[Mc_AN];
1648 Cwan = WhatSpecies[Gc_AN];
1649 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1650 Gh_AN = natn[Gc_AN][h_AN];
1651 Hwan = WhatSpecies[Gh_AN];
1652 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1653 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1654
1655 HisH1[0][0][Mc_AN][h_AN][i][j] = H[0][Mc_AN][h_AN][i][j];
1656 HisH1[0][1][Mc_AN][h_AN][i][j] = H[1][Mc_AN][h_AN][i][j];
1657 HisH1[0][2][Mc_AN][h_AN][i][j] = H[2][Mc_AN][h_AN][i][j];
1658 HisH1[0][3][Mc_AN][h_AN][i][j] = H[3][Mc_AN][h_AN][i][j];
1659
1660 HisH2[0][0][Mc_AN][h_AN][i][j] = iHNL[0][Mc_AN][h_AN][i][j];
1661 HisH2[0][1][Mc_AN][h_AN][i][j] = iHNL[1][Mc_AN][h_AN][i][j];
1662 HisH2[0][2][Mc_AN][h_AN][i][j] = iHNL[2][Mc_AN][h_AN][i][j];
1663 }
1664 }
1665 }
1666 }
1667
1668 }
1669
1670 else {
1671
1672 /* shift the current Hamiltonian */
1673
1674 for (m=dim; 0<m; m--){
1675 for (spin=0; spin<=SpinP_switch; spin++){
1676 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1677 Gc_AN = M2G[Mc_AN];
1678 Cwan = WhatSpecies[Gc_AN];
1679 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1680 Gh_AN = natn[Gc_AN][h_AN];
1681 Hwan = WhatSpecies[Gh_AN];
1682 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1683 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1684 HisH1[m][spin][Mc_AN][h_AN][i][j] = HisH1[m-1][spin][Mc_AN][h_AN][i][j];
1685 }
1686 }
1687 }
1688 }
1689 }
1690 }
1691
1692 /* save the current Hamiltonian */
1693
1694 for (spin=0; spin<=SpinP_switch; spin++){
1695 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1696 Gc_AN = M2G[Mc_AN];
1697 Cwan = WhatSpecies[Gc_AN];
1698 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1699 Gh_AN = natn[Gc_AN][h_AN];
1700 Hwan = WhatSpecies[Gh_AN];
1701 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1702 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1703 HisH1[0][spin][Mc_AN][h_AN][i][j] = H[spin][Mc_AN][h_AN][i][j];
1704 }
1705 }
1706 }
1707 }
1708 }
1709
1710 } /* else */
1711
1712 /****************************************************
1713 freeing of arrays
1714 ****************************************************/
1715
1716 free(coes);
1717
1718 for (i=0; i<List_YOUSO[39]; i++){
1719 free(A[i]);
1720 }
1721 free(A);
1722
1723 for (i=0; i<List_YOUSO[39]; i++){
1724 free(IA[i]);
1725 }
1726 free(IA);
1727 }
1728
1729
1730
1731
1732
Pulay_Mixing_H(int MD_iter,int SCF_iter,int SCF_iter0)1733 void Pulay_Mixing_H(int MD_iter, int SCF_iter, int SCF_iter0 )
1734 {
1735 int Mc_AN,Gc_AN,Cwan,Hwan,h_AN,Gh_AN,i,j,spin;
1736 int dim,m,n,flag_nan,tno;
1737 double my_sum,tmp1,tmp2,alpha,max_diff,d;
1738 double r,r10,r11,r12,r13,r20,r21,r22;
1739 double h,h10,h11,h12,h13,h20,h21,h22;
1740 double **A,**IA,*coes,**metric;
1741 char nanchar[300];
1742
1743 /****************************************************
1744 determination of dimension of the subspace
1745 ****************************************************/
1746
1747 if (SCF_iter<=Num_Mixing_pDM) dim = SCF_iter-1;
1748 else dim = Num_Mixing_pDM;
1749
1750 /****************************************************
1751 allocation of arrays
1752 ****************************************************/
1753
1754 coes = (double*)malloc(sizeof(double)*List_YOUSO[39]);
1755
1756 A = (double**)malloc(sizeof(double*)*List_YOUSO[39]);
1757 for (i=0; i<List_YOUSO[39]; i++){
1758 A[i] = (double*)malloc(sizeof(double)*List_YOUSO[39]);
1759 }
1760
1761 IA = (double**)malloc(sizeof(double*)*List_YOUSO[39]);
1762 for (i=0; i<List_YOUSO[39]; i++){
1763 IA[i] = (double*)malloc(sizeof(double)*List_YOUSO[39]);
1764 }
1765
1766 metric = (double**)malloc(sizeof(double*)*(Matomnum+1));
1767 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
1768 if (Mc_AN==0){
1769 tno = 1;
1770 }
1771 else{
1772 Gc_AN = M2G[Mc_AN];
1773 Cwan = WhatSpecies[Gc_AN];
1774 tno = Spe_Total_NO[Cwan];
1775 }
1776 metric[Mc_AN] = (double*)malloc(sizeof(double)*tno);
1777 }
1778
1779 /****************************************************
1780 determine metric used for calculations of norm
1781 ****************************************************/
1782
1783 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1784 Gc_AN = M2G[Mc_AN];
1785 Cwan = WhatSpecies[Gc_AN];
1786 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1787 d = fabs(HisH1[0][0][Mc_AN][0][i][i]-ChemP);
1788 metric[Mc_AN][i] = 5.0/(d*d+5.0);
1789 }
1790 }
1791
1792 /****************************************************
1793 shift the residual H
1794 ****************************************************/
1795
1796 if (SpinP_switch==3){
1797
1798 /* shift the residual Hamiltonian */
1799
1800 for (m=dim; 0<m; m--){
1801 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1802 Gc_AN = M2G[Mc_AN];
1803 Cwan = WhatSpecies[Gc_AN];
1804 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1805 Gh_AN = natn[Gc_AN][h_AN];
1806 Hwan = WhatSpecies[Gh_AN];
1807 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1808 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1809
1810 ResidualH1[m][0][Mc_AN][h_AN][i][j] = ResidualH1[m-1][0][Mc_AN][h_AN][i][j];
1811 ResidualH1[m][1][Mc_AN][h_AN][i][j] = ResidualH1[m-1][1][Mc_AN][h_AN][i][j];
1812 ResidualH1[m][2][Mc_AN][h_AN][i][j] = ResidualH1[m-1][2][Mc_AN][h_AN][i][j];
1813 ResidualH1[m][3][Mc_AN][h_AN][i][j] = ResidualH1[m-1][3][Mc_AN][h_AN][i][j];
1814
1815 ResidualH2[m][0][Mc_AN][h_AN][i][j] = ResidualH2[m-1][0][Mc_AN][h_AN][i][j];
1816 ResidualH2[m][1][Mc_AN][h_AN][i][j] = ResidualH2[m-1][1][Mc_AN][h_AN][i][j];
1817 ResidualH2[m][2][Mc_AN][h_AN][i][j] = ResidualH2[m-1][2][Mc_AN][h_AN][i][j];
1818 }
1819 }
1820 }
1821 }
1822 }
1823
1824 /* calculate the current residual Hamiltonian */
1825
1826 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1827 Gc_AN = M2G[Mc_AN];
1828 Cwan = WhatSpecies[Gc_AN];
1829 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1830 Gh_AN = natn[Gc_AN][h_AN];
1831 Hwan = WhatSpecies[Gh_AN];
1832 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1833 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1834
1835 ResidualH1[0][0][Mc_AN][h_AN][i][j] = H[0][Mc_AN][h_AN][i][j] - HisH1[0][0][Mc_AN][h_AN][i][j];
1836 ResidualH1[0][1][Mc_AN][h_AN][i][j] = H[1][Mc_AN][h_AN][i][j] - HisH1[0][1][Mc_AN][h_AN][i][j];
1837 ResidualH1[0][2][Mc_AN][h_AN][i][j] = H[2][Mc_AN][h_AN][i][j] - HisH1[0][2][Mc_AN][h_AN][i][j];
1838 ResidualH1[0][3][Mc_AN][h_AN][i][j] = H[3][Mc_AN][h_AN][i][j] - HisH1[0][3][Mc_AN][h_AN][i][j];
1839
1840 ResidualH2[0][0][Mc_AN][h_AN][i][j] = iHNL[0][Mc_AN][h_AN][i][j] - HisH2[0][0][Mc_AN][h_AN][i][j];
1841 ResidualH2[0][1][Mc_AN][h_AN][i][j] = iHNL[1][Mc_AN][h_AN][i][j] - HisH2[0][1][Mc_AN][h_AN][i][j];
1842 ResidualH2[0][2][Mc_AN][h_AN][i][j] = iHNL[2][Mc_AN][h_AN][i][j] - HisH2[0][2][Mc_AN][h_AN][i][j];
1843 }
1844 }
1845 }
1846 }
1847
1848 }
1849
1850 else{
1851
1852 /* shift the residual Hamiltonian */
1853
1854 for (m=dim; 0<m; m--){
1855 for (spin=0; spin<=SpinP_switch; spin++){
1856 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1857 Gc_AN = M2G[Mc_AN];
1858 Cwan = WhatSpecies[Gc_AN];
1859 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1860 Gh_AN = natn[Gc_AN][h_AN];
1861 Hwan = WhatSpecies[Gh_AN];
1862 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1863 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1864 ResidualH1[m][spin][Mc_AN][h_AN][i][j] = ResidualH1[m-1][spin][Mc_AN][h_AN][i][j];
1865 }
1866 }
1867 }
1868 }
1869 }
1870 }
1871
1872 /* calculate the current residual Hamiltonian */
1873
1874 for (spin=0; spin<=SpinP_switch; spin++){
1875 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1876 Gc_AN = M2G[Mc_AN];
1877 Cwan = WhatSpecies[Gc_AN];
1878 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1879 Gh_AN = natn[Gc_AN][h_AN];
1880 Hwan = WhatSpecies[Gh_AN];
1881 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1882 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1883 ResidualH1[0][spin][Mc_AN][h_AN][i][j] = H[spin][Mc_AN][h_AN][i][j] - HisH1[0][spin][Mc_AN][h_AN][i][j];
1884 }
1885 }
1886 }
1887 }
1888 }
1889
1890 } /* else */
1891
1892 /****************************************************
1893 calculation of the residual matrix
1894 ****************************************************/
1895
1896 for (m=0; m<dim; m++){
1897 for (n=0; n<dim; n++){
1898
1899 my_sum = 0.0;
1900
1901 if (SpinP_switch==3){
1902
1903 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1904 Gc_AN = M2G[Mc_AN];
1905 Cwan = WhatSpecies[Gc_AN];
1906 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1907 Gh_AN = natn[Gc_AN][h_AN];
1908 Hwan = WhatSpecies[Gh_AN];
1909 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1910 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1911
1912 tmp1 = ResidualH1[m][0][Mc_AN][h_AN][i][j];
1913 tmp2 = ResidualH1[n][0][Mc_AN][h_AN][i][j];
1914 my_sum += metric[Mc_AN][i]*tmp1*tmp2;
1915
1916 tmp1 = ResidualH1[m][1][Mc_AN][h_AN][i][j];
1917 tmp2 = ResidualH1[n][1][Mc_AN][h_AN][i][j];
1918 my_sum += metric[Mc_AN][i]*tmp1*tmp2;
1919
1920 tmp1 = ResidualH1[m][2][Mc_AN][h_AN][i][j];
1921 tmp2 = ResidualH1[n][2][Mc_AN][h_AN][i][j];
1922 my_sum += metric[Mc_AN][i]*tmp1*tmp2;
1923
1924 tmp1 = ResidualH1[m][3][Mc_AN][h_AN][i][j];
1925 tmp2 = ResidualH1[n][3][Mc_AN][h_AN][i][j];
1926 my_sum += metric[Mc_AN][i]*tmp1*tmp2;
1927
1928 tmp1 = ResidualH2[m][0][Mc_AN][h_AN][i][j];
1929 tmp2 = ResidualH2[n][0][Mc_AN][h_AN][i][j];
1930 my_sum += metric[Mc_AN][i]*tmp1*tmp2;
1931
1932 tmp1 = ResidualH2[m][1][Mc_AN][h_AN][i][j];
1933 tmp2 = ResidualH2[n][1][Mc_AN][h_AN][i][j];
1934 my_sum += metric[Mc_AN][i]*tmp1*tmp2;
1935
1936 tmp1 = ResidualH2[m][2][Mc_AN][h_AN][i][j];
1937 tmp2 = ResidualH2[n][2][Mc_AN][h_AN][i][j];
1938 my_sum += metric[Mc_AN][i]*tmp1*tmp2;
1939 }
1940 }
1941 }
1942 }
1943
1944 } /* if (SpinP_switch==3 */
1945
1946 else{
1947
1948 for (spin=0; spin<=SpinP_switch; spin++){
1949 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1950 Gc_AN = M2G[Mc_AN];
1951 Cwan = WhatSpecies[Gc_AN];
1952 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1953 Gh_AN = natn[Gc_AN][h_AN];
1954 Hwan = WhatSpecies[Gh_AN];
1955 for (i=0; i<Spe_Total_NO[Cwan]; i++){
1956 for (j=0; j<Spe_Total_NO[Hwan]; j++){
1957 tmp1 = ResidualH1[m][spin][Mc_AN][h_AN][i][j];
1958 tmp2 = ResidualH1[n][spin][Mc_AN][h_AN][i][j];
1959 my_sum += metric[Mc_AN][i]*tmp1*tmp2;
1960 }
1961 }
1962 }
1963 }
1964 }
1965
1966 } /* else */
1967
1968 MPI_Allreduce(&my_sum, &A[m][n], 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
1969 A[n][m] = A[m][n];
1970
1971 } /* n */
1972 } /* m */
1973
1974 NormRD[0] = A[0][0];
1975
1976 for (m=1; m<=dim; m++){
1977 A[m-1][dim] = -1.0;
1978 A[dim][m-1] = -1.0;
1979 }
1980 A[dim][dim] = 0.0;
1981
1982 Inverse(dim,A,IA);
1983
1984 for (m=1; m<=dim; m++){
1985 coes[m] = -IA[m-1][dim];
1986 }
1987
1988 /****************************************************
1989 check "nan", "NaN", "inf" or "Inf"
1990 ****************************************************/
1991
1992 flag_nan = 0;
1993 for (m=1; m<=dim; m++){
1994
1995 sprintf(nanchar,"%8.4f",coes[m]);
1996 if ( strstr(nanchar,"nan")!=NULL || strstr(nanchar,"NaN")!=NULL
1997 || strstr(nanchar,"inf")!=NULL || strstr(nanchar,"Inf")!=NULL){
1998
1999 flag_nan = 1;
2000 }
2001 }
2002
2003 if (flag_nan==1){
2004 for (m=1; m<=dim; m++){
2005 coes[m] = 0.0;
2006 }
2007 coes[1] = 0.05;
2008 coes[2] = 0.95;
2009 }
2010
2011 /****************************************************
2012 calculation of optimum residual Hamiltonian
2013 ****************************************************/
2014
2015 if (SpinP_switch==3){
2016
2017 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2018 Gc_AN = M2G[Mc_AN];
2019 Cwan = WhatSpecies[Gc_AN];
2020 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2021 Gh_AN = natn[Gc_AN][h_AN];
2022 Hwan = WhatSpecies[Gh_AN];
2023 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2024 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2025
2026 r10 = 0.0;
2027 r11 = 0.0;
2028 r12 = 0.0;
2029 r13 = 0.0;
2030
2031 r20 = 0.0;
2032 r21 = 0.0;
2033 r22 = 0.0;
2034
2035 for (m=0; m<dim; m++){
2036
2037 r10 += ResidualH1[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
2038 r11 += ResidualH1[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
2039 r12 += ResidualH1[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
2040 r13 += ResidualH1[m][3][Mc_AN][h_AN][i][j]*coes[m+1];
2041
2042 r20 += ResidualH2[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
2043 r21 += ResidualH2[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
2044 r22 += ResidualH2[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
2045 }
2046
2047 /* optimum Residual H is stored in ResidualH1[dim] and ResidualH2[dim] */
2048
2049 ResidualH1[dim][0][Mc_AN][h_AN][i][j] = r10;
2050 ResidualH1[dim][1][Mc_AN][h_AN][i][j] = r11;
2051 ResidualH1[dim][2][Mc_AN][h_AN][i][j] = r12;
2052 ResidualH1[dim][3][Mc_AN][h_AN][i][j] = r13;
2053
2054 ResidualH2[dim][0][Mc_AN][h_AN][i][j] = r20;
2055 ResidualH2[dim][1][Mc_AN][h_AN][i][j] = r21;
2056 ResidualH2[dim][2][Mc_AN][h_AN][i][j] = r22;
2057
2058 }
2059 }
2060 }
2061 }
2062
2063 }
2064
2065 else{
2066
2067 for (spin=0; spin<=SpinP_switch; spin++){
2068 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2069 Gc_AN = M2G[Mc_AN];
2070 Cwan = WhatSpecies[Gc_AN];
2071 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2072 Gh_AN = natn[Gc_AN][h_AN];
2073 Hwan = WhatSpecies[Gh_AN];
2074 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2075 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2076
2077 r = 0.0;
2078 for (m=0; m<dim; m++){
2079 r += ResidualH1[m][spin][Mc_AN][h_AN][i][j]*coes[m+1];
2080 }
2081
2082 /* optimum Residual H is stored in ResidualH1[dim] */
2083
2084 ResidualH1[dim][spin][Mc_AN][h_AN][i][j] = r;
2085
2086 }
2087 }
2088 }
2089 }
2090 }
2091 }
2092
2093 /****************************************************
2094 mixing of Hamiltonian
2095 ****************************************************/
2096
2097 if (1.0e-1<=NormRD[0])
2098 alpha = 0.5;
2099 else if (1.0e-2<=NormRD[0] && NormRD[0]<1.0e-1)
2100 alpha = 0.6;
2101 else if (1.0e-3<=NormRD[0] && NormRD[0]<1.0e-2)
2102 alpha = 0.7;
2103 else if (1.0e-4<=NormRD[0] && NormRD[0]<1.0e-3)
2104 alpha = 0.8;
2105 else
2106 alpha = 1.0;
2107
2108 if (SpinP_switch==3){
2109
2110 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2111 Gc_AN = M2G[Mc_AN];
2112 Cwan = WhatSpecies[Gc_AN];
2113 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2114 Gh_AN = natn[Gc_AN][h_AN];
2115 Hwan = WhatSpecies[Gh_AN];
2116 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2117 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2118
2119 h10 = 0.0;
2120 h11 = 0.0;
2121 h12 = 0.0;
2122 h13 = 0.0;
2123
2124 h20 = 0.0;
2125 h21 = 0.0;
2126 h22 = 0.0;
2127
2128 for (m=0; m<dim; m++){
2129
2130 h10 += HisH1[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
2131 h11 += HisH1[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
2132 h12 += HisH1[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
2133 h13 += HisH1[m][3][Mc_AN][h_AN][i][j]*coes[m+1];
2134
2135 h20 += HisH2[m][0][Mc_AN][h_AN][i][j]*coes[m+1];
2136 h21 += HisH2[m][1][Mc_AN][h_AN][i][j]*coes[m+1];
2137 h22 += HisH2[m][2][Mc_AN][h_AN][i][j]*coes[m+1];
2138 }
2139
2140 H[0][Mc_AN][h_AN][i][j] = h10 + alpha*ResidualH1[dim][0][Mc_AN][h_AN][i][j];
2141 H[1][Mc_AN][h_AN][i][j] = h11 + alpha*ResidualH1[dim][1][Mc_AN][h_AN][i][j];
2142 H[2][Mc_AN][h_AN][i][j] = h12 + alpha*ResidualH1[dim][2][Mc_AN][h_AN][i][j];
2143 H[3][Mc_AN][h_AN][i][j] = h13 + alpha*ResidualH1[dim][3][Mc_AN][h_AN][i][j];
2144
2145 iHNL[0][Mc_AN][h_AN][i][j] = h20 + alpha*ResidualH2[dim][0][Mc_AN][h_AN][i][j];
2146 iHNL[1][Mc_AN][h_AN][i][j] = h21 + alpha*ResidualH2[dim][1][Mc_AN][h_AN][i][j];
2147 iHNL[2][Mc_AN][h_AN][i][j] = h22 + alpha*ResidualH2[dim][2][Mc_AN][h_AN][i][j];
2148
2149 }
2150 }
2151 }
2152 }
2153
2154 }
2155
2156 else{
2157
2158 for (spin=0; spin<=SpinP_switch; spin++){
2159 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2160 Gc_AN = M2G[Mc_AN];
2161 Cwan = WhatSpecies[Gc_AN];
2162 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2163 Gh_AN = natn[Gc_AN][h_AN];
2164 Hwan = WhatSpecies[Gh_AN];
2165 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2166 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2167
2168 r = 0.0;
2169 h = 0.0;
2170
2171 for (m=0; m<dim; m++){
2172 r += ResidualH1[m][spin][Mc_AN][h_AN][i][j]*coes[m+1];
2173 h += HisH1[m][spin][Mc_AN][h_AN][i][j]*coes[m+1];
2174 }
2175
2176 H[spin][Mc_AN][h_AN][i][j] = h + alpha*r;
2177 }
2178 }
2179 }
2180 }
2181 }
2182 }
2183
2184 /****************************************************
2185 shifting of Hamiltonian
2186 ****************************************************/
2187
2188 if (SpinP_switch==3){
2189
2190 /* shift the current Hamiltonian */
2191
2192 for (m=dim; 0<m; m--){
2193 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2194 Gc_AN = M2G[Mc_AN];
2195 Cwan = WhatSpecies[Gc_AN];
2196 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2197 Gh_AN = natn[Gc_AN][h_AN];
2198 Hwan = WhatSpecies[Gh_AN];
2199 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2200 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2201
2202 HisH1[m][0][Mc_AN][h_AN][i][j] = HisH1[m-1][0][Mc_AN][h_AN][i][j];
2203 HisH1[m][1][Mc_AN][h_AN][i][j] = HisH1[m-1][1][Mc_AN][h_AN][i][j];
2204 HisH1[m][2][Mc_AN][h_AN][i][j] = HisH1[m-1][2][Mc_AN][h_AN][i][j];
2205 HisH1[m][3][Mc_AN][h_AN][i][j] = HisH1[m-1][3][Mc_AN][h_AN][i][j];
2206
2207 HisH2[m][0][Mc_AN][h_AN][i][j] = HisH2[m-1][0][Mc_AN][h_AN][i][j];
2208 HisH2[m][1][Mc_AN][h_AN][i][j] = HisH2[m-1][1][Mc_AN][h_AN][i][j];
2209 HisH2[m][2][Mc_AN][h_AN][i][j] = HisH2[m-1][2][Mc_AN][h_AN][i][j];
2210 }
2211 }
2212 }
2213 }
2214 }
2215
2216 /* save the current Hamiltonian */
2217
2218 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2219 Gc_AN = M2G[Mc_AN];
2220 Cwan = WhatSpecies[Gc_AN];
2221 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2222 Gh_AN = natn[Gc_AN][h_AN];
2223 Hwan = WhatSpecies[Gh_AN];
2224 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2225 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2226
2227 HisH1[0][0][Mc_AN][h_AN][i][j] = H[0][Mc_AN][h_AN][i][j];
2228 HisH1[0][1][Mc_AN][h_AN][i][j] = H[1][Mc_AN][h_AN][i][j];
2229 HisH1[0][2][Mc_AN][h_AN][i][j] = H[2][Mc_AN][h_AN][i][j];
2230 HisH1[0][3][Mc_AN][h_AN][i][j] = H[3][Mc_AN][h_AN][i][j];
2231
2232 HisH2[0][0][Mc_AN][h_AN][i][j] = iHNL[0][Mc_AN][h_AN][i][j];
2233 HisH2[0][1][Mc_AN][h_AN][i][j] = iHNL[1][Mc_AN][h_AN][i][j];
2234 HisH2[0][2][Mc_AN][h_AN][i][j] = iHNL[2][Mc_AN][h_AN][i][j];
2235 }
2236 }
2237 }
2238 }
2239
2240 }
2241
2242 else {
2243
2244 /* shift the current Hamiltonian */
2245
2246 for (m=dim; 0<m; m--){
2247 for (spin=0; spin<=SpinP_switch; spin++){
2248 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2249 Gc_AN = M2G[Mc_AN];
2250 Cwan = WhatSpecies[Gc_AN];
2251 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2252 Gh_AN = natn[Gc_AN][h_AN];
2253 Hwan = WhatSpecies[Gh_AN];
2254 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2255 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2256 HisH1[m][spin][Mc_AN][h_AN][i][j] = HisH1[m-1][spin][Mc_AN][h_AN][i][j];
2257 }
2258 }
2259 }
2260 }
2261 }
2262 }
2263
2264 /* save the current Hamiltonian */
2265
2266 for (spin=0; spin<=SpinP_switch; spin++){
2267 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2268 Gc_AN = M2G[Mc_AN];
2269 Cwan = WhatSpecies[Gc_AN];
2270 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2271 Gh_AN = natn[Gc_AN][h_AN];
2272 Hwan = WhatSpecies[Gh_AN];
2273 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2274 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2275 HisH1[0][spin][Mc_AN][h_AN][i][j] = H[spin][Mc_AN][h_AN][i][j];
2276 }
2277 }
2278 }
2279 }
2280 }
2281
2282 } /* else */
2283
2284 /****************************************************
2285 freeing of arrays
2286 ****************************************************/
2287
2288 free(coes);
2289
2290 for (i=0; i<List_YOUSO[39]; i++){
2291 free(A[i]);
2292 }
2293 free(A);
2294
2295 for (i=0; i<List_YOUSO[39]; i++){
2296 free(IA[i]);
2297 }
2298 free(IA);
2299
2300 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
2301 if (Mc_AN==0){
2302 tno = 1;
2303 }
2304 else{
2305 Gc_AN = M2G[Mc_AN];
2306 Cwan = WhatSpecies[Gc_AN];
2307 tno = Spe_Total_NO[Cwan];
2308 }
2309 free(metric[Mc_AN]);
2310 }
2311 free(metric);
2312
2313 }
2314
2315
2316
Simple_Mixing_H(int MD_iter,int SCF_iter,int SCF_iter0)2317 void Simple_Mixing_H(int MD_iter, int SCF_iter, int SCF_iter0 )
2318 {
2319 int m,Mc_AN,Gc_AN,Cwan,spin,dim;
2320 int i,j,h_AN,Gh_AN,Hwan,ian,jan,n;
2321 double w1,w2,My_Norm,Norm;
2322 double d0,d1,d2,d3,d,tmp0;
2323 double Mix_wgt,Min_Weight,Max_Weight;
2324 int numprocs,myid,ID;
2325
2326 /* MPI */
2327 MPI_Comm_size(mpi_comm_level1,&numprocs);
2328 MPI_Comm_rank(mpi_comm_level1,&myid);
2329
2330 /* start... */
2331
2332 Min_Weight = Min_Mixing_weight;
2333 if (SCF_RENZOKU==-1){
2334 Max_Weight = Max_Mixing_weight;
2335 Max_Mixing_weight2 = Max_Mixing_weight;
2336 }
2337 else if (SCF_RENZOKU==1000){ /* past 3 */
2338 Max_Mixing_weight2 = 2.0*Max_Mixing_weight2;
2339 if (0.7<Max_Mixing_weight2) Max_Mixing_weight2 = 0.7;
2340 Max_Weight = Max_Mixing_weight2;
2341 SCF_RENZOKU = 0;
2342 }
2343 else{
2344 Max_Weight = Max_Mixing_weight2;
2345 }
2346
2347 /* determination of dim */
2348
2349 if (SCF_iter<Num_Mixing_pDM) dim = SCF_iter;
2350 else dim = Num_Mixing_pDM;
2351
2352 /****************************************************
2353 shift the residual H
2354 ****************************************************/
2355
2356 if (Mixing_switch!=1 && Mixing_switch!=6){
2357
2358 if (SpinP_switch==3){
2359
2360 /* shift the residual Hamiltonian */
2361
2362 for (m=(dim-1); 0<m; m--){
2363 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2364 Gc_AN = M2G[Mc_AN];
2365 Cwan = WhatSpecies[Gc_AN];
2366 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2367 Gh_AN = natn[Gc_AN][h_AN];
2368 Hwan = WhatSpecies[Gh_AN];
2369 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2370 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2371
2372 ResidualH1[m][0][Mc_AN][h_AN][i][j] = ResidualH1[m-1][0][Mc_AN][h_AN][i][j];
2373 ResidualH1[m][1][Mc_AN][h_AN][i][j] = ResidualH1[m-1][1][Mc_AN][h_AN][i][j];
2374 ResidualH1[m][2][Mc_AN][h_AN][i][j] = ResidualH1[m-1][2][Mc_AN][h_AN][i][j];
2375 ResidualH1[m][3][Mc_AN][h_AN][i][j] = ResidualH1[m-1][3][Mc_AN][h_AN][i][j];
2376
2377 ResidualH2[m][0][Mc_AN][h_AN][i][j] = ResidualH2[m-1][0][Mc_AN][h_AN][i][j];
2378 ResidualH2[m][1][Mc_AN][h_AN][i][j] = ResidualH2[m-1][1][Mc_AN][h_AN][i][j];
2379 ResidualH2[m][2][Mc_AN][h_AN][i][j] = ResidualH2[m-1][2][Mc_AN][h_AN][i][j];
2380 }
2381 }
2382 }
2383 }
2384 }
2385
2386 /* calculate the current residual Hamiltonian */
2387
2388 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2389 Gc_AN = M2G[Mc_AN];
2390 Cwan = WhatSpecies[Gc_AN];
2391 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2392 Gh_AN = natn[Gc_AN][h_AN];
2393 Hwan = WhatSpecies[Gh_AN];
2394 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2395 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2396
2397 ResidualH1[0][0][Mc_AN][h_AN][i][j] = H[0][Mc_AN][h_AN][i][j] - HisH1[0][0][Mc_AN][h_AN][i][j];
2398 ResidualH1[0][1][Mc_AN][h_AN][i][j] = H[1][Mc_AN][h_AN][i][j] - HisH1[0][1][Mc_AN][h_AN][i][j];
2399 ResidualH1[0][2][Mc_AN][h_AN][i][j] = H[2][Mc_AN][h_AN][i][j] - HisH1[0][2][Mc_AN][h_AN][i][j];
2400 ResidualH1[0][3][Mc_AN][h_AN][i][j] = H[3][Mc_AN][h_AN][i][j] - HisH1[0][3][Mc_AN][h_AN][i][j];
2401
2402 ResidualH2[0][0][Mc_AN][h_AN][i][j] = iHNL[0][Mc_AN][h_AN][i][j] - HisH2[0][0][Mc_AN][h_AN][i][j];
2403 ResidualH2[0][1][Mc_AN][h_AN][i][j] = iHNL[1][Mc_AN][h_AN][i][j] - HisH2[0][1][Mc_AN][h_AN][i][j];
2404 ResidualH2[0][2][Mc_AN][h_AN][i][j] = iHNL[2][Mc_AN][h_AN][i][j] - HisH2[0][2][Mc_AN][h_AN][i][j];
2405 }
2406 }
2407 }
2408 }
2409
2410 }
2411
2412 else{
2413
2414 /* shift the residual Hamiltonian */
2415
2416 for (m=(dim-1); 0<m; m--){
2417 for (spin=0; spin<=SpinP_switch; spin++){
2418 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2419 Gc_AN = M2G[Mc_AN];
2420 Cwan = WhatSpecies[Gc_AN];
2421 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2422 Gh_AN = natn[Gc_AN][h_AN];
2423 Hwan = WhatSpecies[Gh_AN];
2424 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2425 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2426 ResidualH1[m][spin][Mc_AN][h_AN][i][j] = ResidualH1[m-1][spin][Mc_AN][h_AN][i][j];
2427 }
2428 }
2429 }
2430 }
2431 }
2432 }
2433
2434 /* calculate the current residual Hamiltonian */
2435
2436 for (spin=0; spin<=SpinP_switch; spin++){
2437 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2438 Gc_AN = M2G[Mc_AN];
2439 Cwan = WhatSpecies[Gc_AN];
2440 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2441 Gh_AN = natn[Gc_AN][h_AN];
2442 Hwan = WhatSpecies[Gh_AN];
2443 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2444 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2445 ResidualH1[0][spin][Mc_AN][h_AN][i][j] = H[spin][Mc_AN][h_AN][i][j] - HisH1[0][spin][Mc_AN][h_AN][i][j];
2446 }
2447 }
2448 }
2449 }
2450 }
2451
2452 } /* else */
2453
2454 } /* if (Mixing_switch!=1 && Mixing_switch!=6) */
2455
2456 /****************************************************
2457 calculation of NormRH
2458 ****************************************************/
2459
2460 My_Norm = 0.0;
2461
2462 if (SpinP_switch==3){
2463
2464 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2465 Gc_AN = M2G[Mc_AN];
2466 Cwan = WhatSpecies[Gc_AN];
2467 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2468 Gh_AN = natn[Gc_AN][h_AN];
2469 Hwan = WhatSpecies[Gh_AN];
2470 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2471 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2472
2473 d0 = HisH1[0][0][Mc_AN][h_AN][i][j] - H[0][Mc_AN][h_AN][i][j];
2474 d1 = HisH1[0][1][Mc_AN][h_AN][i][j] - H[1][Mc_AN][h_AN][i][j];
2475 d2 = HisH1[0][2][Mc_AN][h_AN][i][j] - H[2][Mc_AN][h_AN][i][j];
2476 d3 = HisH1[0][3][Mc_AN][h_AN][i][j] - H[3][Mc_AN][h_AN][i][j];
2477 My_Norm += d0*d0 + d1*d1 + d2*d2 + d3*d3;
2478
2479 d0 = HisH2[0][0][Mc_AN][h_AN][i][j] - iHNL[0][Mc_AN][h_AN][i][j];
2480 d1 = HisH2[0][1][Mc_AN][h_AN][i][j] - iHNL[1][Mc_AN][h_AN][i][j];
2481 d2 = HisH2[0][2][Mc_AN][h_AN][i][j] - iHNL[2][Mc_AN][h_AN][i][j];
2482 My_Norm += d0*d0 + d1*d1 + d2*d2;
2483
2484 }
2485 }
2486 }
2487 }
2488
2489 }
2490
2491 else{
2492
2493 for (spin=0; spin<=SpinP_switch; spin++){
2494 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2495 Gc_AN = M2G[Mc_AN];
2496 Cwan = WhatSpecies[Gc_AN];
2497 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2498 Gh_AN = natn[Gc_AN][h_AN];
2499 Hwan = WhatSpecies[Gh_AN];
2500 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2501 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2502 d = HisH1[0][spin][Mc_AN][h_AN][i][j] - H[spin][Mc_AN][h_AN][i][j];
2503 My_Norm += d*d;
2504 }
2505 }
2506 }
2507 }
2508 }
2509 }
2510
2511 /****************************************************
2512 MPI:
2513
2514 My_Norm
2515 ****************************************************/
2516
2517 MPI_Allreduce(&My_Norm, &Norm, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
2518
2519 /****************************************************
2520 find an optimum mixing weight
2521 ****************************************************/
2522
2523 for (i=4; 1<=i; i--){
2524 NormRD[i] = NormRD[i-1];
2525 History_Uele[i] = History_Uele[i-1];
2526 }
2527 NormRD[0] = Norm;
2528 History_Uele[0] = Uele;
2529
2530 if (1<SCF_iter){
2531
2532 if ( (int)sgn(History_Uele[0]-History_Uele[1])
2533 ==(int)sgn(History_Uele[1]-History_Uele[2])
2534 && NormRD[0]<NormRD[1]){
2535
2536 /* tmp0 = 1.6*Mixing_weight; */
2537
2538 tmp0 = NormRD[1]/(largest(NormRD[1]-NormRD[0],10e-10))*Mixing_weight;
2539
2540 if (tmp0<Max_Weight){
2541 if (Min_Weight<tmp0){
2542 Mixing_weight = tmp0;
2543 }
2544 else{
2545 Mixing_weight = Min_Weight;
2546 }
2547 }
2548 else{
2549 Mixing_weight = Max_Weight;
2550 SCF_RENZOKU++;
2551 }
2552 }
2553
2554 else if ( (int)sgn(History_Uele[0]-History_Uele[1])
2555 ==(int)sgn(History_Uele[1]-History_Uele[2])
2556 && NormRD[1]<NormRD[0]){
2557
2558 tmp0 = NormRD[1]/(largest(NormRD[1]+NormRD[0],10e-10))*Mixing_weight;
2559
2560 /* tmp0 = Mixing_weight/1.6; */
2561
2562 if (tmp0<Max_Weight){
2563 if (Min_Weight<tmp0)
2564 Mixing_weight = tmp0;
2565 else
2566 Mixing_weight = Min_Weight;
2567 }
2568 else{
2569 Mixing_weight = Max_Weight;
2570 }
2571
2572 SCF_RENZOKU = -1;
2573 }
2574
2575 else if ( (int)sgn(History_Uele[0]-History_Uele[1])
2576 !=(int)sgn(History_Uele[1]-History_Uele[2])
2577 && NormRD[0]<NormRD[1]){
2578
2579 /* tmp0 = Mixing_weight*1.2; */
2580
2581 tmp0 = NormRD[1]/(largest(NormRD[1]-NormRD[0],10e-10))*Mixing_weight;
2582
2583 if (tmp0<Max_Weight){
2584 if (Min_Weight<tmp0)
2585 Mixing_weight = tmp0;
2586 else
2587 Mixing_weight = Min_Weight;
2588 }
2589 else{
2590 Mixing_weight = Max_Weight;
2591 SCF_RENZOKU++;
2592 }
2593 }
2594
2595 else if ( (int)sgn(History_Uele[0]-History_Uele[1])
2596 !=(int)sgn(History_Uele[1]-History_Uele[2])
2597 && NormRD[1]<NormRD[0]){
2598
2599 /* tmp0 = Mixing_weight/2.0; */
2600
2601 tmp0 = NormRD[1]/(largest(NormRD[1]+NormRD[0],10e-10))*Mixing_weight;
2602
2603 if (tmp0<Max_Weight){
2604 if (Min_Weight<tmp0)
2605 Mixing_weight = tmp0;
2606 else
2607 Mixing_weight = Min_Weight;
2608 }
2609 else
2610 Mixing_weight = Max_Weight;
2611
2612 SCF_RENZOKU = -1;
2613 }
2614 }
2615
2616 Mix_wgt = Mixing_weight;
2617
2618 /****************************************************
2619 finding a proper mixing weight
2620 ****************************************************/
2621
2622 if (SCF_iter==1){
2623 w1 = 1.0;
2624 w2 = 1.0 - w1;
2625 }
2626 else{
2627 w1 = Mix_wgt;
2628 w2 = 1.0 - w1;
2629 }
2630
2631 /****************************************************
2632 performing the simple mixing for Hamiltonian
2633 ****************************************************/
2634
2635 if (SpinP_switch==3){
2636
2637 /* shift the current Hamiltonian */
2638
2639 for (m=dim; 0<m; m--){
2640 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2641 Gc_AN = M2G[Mc_AN];
2642 Cwan = WhatSpecies[Gc_AN];
2643 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2644 Gh_AN = natn[Gc_AN][h_AN];
2645 Hwan = WhatSpecies[Gh_AN];
2646 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2647 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2648
2649 HisH1[m][0][Mc_AN][h_AN][i][j] = HisH1[m-1][0][Mc_AN][h_AN][i][j];
2650 HisH1[m][1][Mc_AN][h_AN][i][j] = HisH1[m-1][1][Mc_AN][h_AN][i][j];
2651 HisH1[m][2][Mc_AN][h_AN][i][j] = HisH1[m-1][2][Mc_AN][h_AN][i][j];
2652 HisH1[m][3][Mc_AN][h_AN][i][j] = HisH1[m-1][3][Mc_AN][h_AN][i][j];
2653
2654 HisH2[m][0][Mc_AN][h_AN][i][j] = HisH2[m-1][0][Mc_AN][h_AN][i][j];
2655 HisH2[m][1][Mc_AN][h_AN][i][j] = HisH2[m-1][1][Mc_AN][h_AN][i][j];
2656 HisH2[m][2][Mc_AN][h_AN][i][j] = HisH2[m-1][2][Mc_AN][h_AN][i][j];
2657 }
2658 }
2659 }
2660 }
2661 }
2662
2663 /* mix the current Hamiltonian and the last one */
2664
2665 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2666 Gc_AN = M2G[Mc_AN];
2667 Cwan = WhatSpecies[Gc_AN];
2668 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2669 Gh_AN = natn[Gc_AN][h_AN];
2670 Hwan = WhatSpecies[Gh_AN];
2671 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2672 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2673
2674 H[0][Mc_AN][h_AN][i][j] = w2*HisH1[1][0][Mc_AN][h_AN][i][j] + w1*H[0][Mc_AN][h_AN][i][j];
2675 H[1][Mc_AN][h_AN][i][j] = w2*HisH1[1][1][Mc_AN][h_AN][i][j] + w1*H[1][Mc_AN][h_AN][i][j];
2676 H[2][Mc_AN][h_AN][i][j] = w2*HisH1[1][2][Mc_AN][h_AN][i][j] + w1*H[2][Mc_AN][h_AN][i][j];
2677 H[3][Mc_AN][h_AN][i][j] = w2*HisH1[1][3][Mc_AN][h_AN][i][j] + w1*H[3][Mc_AN][h_AN][i][j];
2678
2679 HisH1[0][0][Mc_AN][h_AN][i][j] = H[0][Mc_AN][h_AN][i][j];
2680 HisH1[0][1][Mc_AN][h_AN][i][j] = H[1][Mc_AN][h_AN][i][j];
2681 HisH1[0][2][Mc_AN][h_AN][i][j] = H[2][Mc_AN][h_AN][i][j];
2682 HisH1[0][3][Mc_AN][h_AN][i][j] = H[3][Mc_AN][h_AN][i][j];
2683
2684 iHNL[0][Mc_AN][h_AN][i][j] = w2*HisH2[1][0][Mc_AN][h_AN][i][j] + w1*iHNL[0][Mc_AN][h_AN][i][j];
2685 iHNL[1][Mc_AN][h_AN][i][j] = w2*HisH2[1][1][Mc_AN][h_AN][i][j] + w1*iHNL[1][Mc_AN][h_AN][i][j];
2686 iHNL[2][Mc_AN][h_AN][i][j] = w2*HisH2[1][2][Mc_AN][h_AN][i][j] + w1*iHNL[2][Mc_AN][h_AN][i][j];
2687
2688 HisH2[0][0][Mc_AN][h_AN][i][j] = iHNL[0][Mc_AN][h_AN][i][j];
2689 HisH2[0][1][Mc_AN][h_AN][i][j] = iHNL[1][Mc_AN][h_AN][i][j];
2690 HisH2[0][2][Mc_AN][h_AN][i][j] = iHNL[2][Mc_AN][h_AN][i][j];
2691 }
2692 }
2693 }
2694 }
2695
2696 }
2697
2698 else {
2699
2700 /* shift the current Hamiltonian */
2701
2702 for (m=dim; 0<m; m--){
2703 for (spin=0; spin<=SpinP_switch; spin++){
2704 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2705 Gc_AN = M2G[Mc_AN];
2706 Cwan = WhatSpecies[Gc_AN];
2707 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2708 Gh_AN = natn[Gc_AN][h_AN];
2709 Hwan = WhatSpecies[Gh_AN];
2710 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2711 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2712 HisH1[m][spin][Mc_AN][h_AN][i][j] = HisH1[m-1][spin][Mc_AN][h_AN][i][j];
2713 }
2714 }
2715 }
2716 }
2717 }
2718 }
2719
2720 /* mix the current Hamiltonian and the last one */
2721
2722 for (spin=0; spin<=SpinP_switch; spin++){
2723 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2724 Gc_AN = M2G[Mc_AN];
2725 Cwan = WhatSpecies[Gc_AN];
2726 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2727 Gh_AN = natn[Gc_AN][h_AN];
2728 Hwan = WhatSpecies[Gh_AN];
2729 for (i=0; i<Spe_Total_NO[Cwan]; i++){
2730 for (j=0; j<Spe_Total_NO[Hwan]; j++){
2731 H[spin][Mc_AN][h_AN][i][j] = w2*HisH1[1][spin][Mc_AN][h_AN][i][j] + w1*H[spin][Mc_AN][h_AN][i][j];
2732 HisH1[0][spin][Mc_AN][h_AN][i][j] = H[spin][Mc_AN][h_AN][i][j];
2733 }
2734 }
2735 }
2736 }
2737 }
2738
2739 } /* else */
2740
2741 /* In case of RMM-DIIS */
2742
2743 if (Mixing_switch==1 || Mixing_switch==6){
2744
2745 for (spin=0; spin<=SpinP_switch; spin++){
2746
2747 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2748 Gc_AN = M2G[Mc_AN];
2749 ian = Spe_Total_CNO[WhatSpecies[Gc_AN]];
2750 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2751 Gh_AN = natn[Gc_AN][h_AN];
2752 jan = Spe_Total_CNO[WhatSpecies[Gh_AN]];
2753
2754 for (m=0; m<ian; m++){
2755 for (n=0; n<jan; n++){
2756
2757 ResidualDM[2][spin][Mc_AN][h_AN][m][n] = DM[0][spin][Mc_AN][h_AN][m][n]
2758 -DM[1][spin][Mc_AN][h_AN][m][n];
2759
2760 DM[2][spin][Mc_AN][h_AN][m][n] = DM[1][spin][Mc_AN][h_AN][m][n];
2761 DM[1][spin][Mc_AN][h_AN][m][n] = DM[0][spin][Mc_AN][h_AN][m][n];
2762 }
2763 }
2764
2765 if ( (SpinP_switch==3 && ( SO_switch==1 || Hub_U_switch==1 || 1<=Constraint_NCS_switch
2766 || Zeeman_NCS_switch==1 || Zeeman_NCO_switch==1 )) && spin<=1 ){
2767
2768 for (m=0; m<ian; m++){
2769 for (n=0; n<jan; n++){
2770
2771 iResidualDM[2][spin][Mc_AN][h_AN][m][n] = iDM[0][spin][Mc_AN][h_AN][m][n]
2772 -iDM[1][spin][Mc_AN][h_AN][m][n];
2773
2774 iDM[2][spin][Mc_AN][h_AN][m][n] = iDM[1][spin][Mc_AN][h_AN][m][n];
2775 iDM[1][spin][Mc_AN][h_AN][m][n] = iDM[0][spin][Mc_AN][h_AN][m][n];
2776 }
2777 }
2778 }
2779
2780 }
2781 }
2782
2783 } /* spin */
2784
2785 } /* if (Mixing_switch==1 || Mixing_switch==6) */
2786
2787 }
2788
2789
2790
2791
2792
2793
Inverse(int n,double ** a,double ** ia)2794 void Inverse(int n, double **a, double **ia)
2795 {
2796 /****************************************************
2797 LU decomposition
2798 0 to n
2799 NOTE:
2800 This routine does consider the reduction of rank
2801 ****************************************************/
2802
2803 int i,j,k;
2804 double w;
2805 double *x,*y;
2806 double **da;
2807
2808 /***************************************************
2809 allocation of arrays:
2810
2811 x[List_YOUSO[39]]
2812 y[List_YOUSO[39]]
2813 da[List_YOUSO[39]][List_YOUSO[39]]
2814 ***************************************************/
2815
2816 x = (double*)malloc(sizeof(double)*List_YOUSO[39]);
2817 y = (double*)malloc(sizeof(double)*List_YOUSO[39]);
2818
2819 da = (double**)malloc(sizeof(double*)*List_YOUSO[39]);
2820 for (i=0; i<List_YOUSO[39]; i++){
2821 da[i] = (double*)malloc(sizeof(double)*List_YOUSO[39]);
2822 }
2823
2824 /* start calc. */
2825
2826 if (n==-1){
2827 for (i=0; i<List_YOUSO[39]; i++){
2828 for (j=0; j<List_YOUSO[39]; j++){
2829 a[i][j] = 0.0;
2830 }
2831 }
2832 }
2833 else{
2834 for (i=0; i<=n; i++){
2835 for (j=0; j<=n; j++){
2836 da[i][j] = a[i][j];
2837 }
2838 }
2839
2840 /****************************************************
2841 LU factorization
2842 ****************************************************/
2843
2844 for (k=0; k<=n-1; k++){
2845 w = 1.0/a[k][k];
2846 for (i=k+1; i<=n; i++){
2847 a[i][k] = w*a[i][k];
2848 for (j=k+1; j<=n; j++){
2849 a[i][j] = a[i][j] - a[i][k]*a[k][j];
2850 }
2851 }
2852 }
2853 for (k=0; k<=n; k++){
2854
2855 /****************************************************
2856 Ly = b
2857 ****************************************************/
2858
2859 for (i=0; i<=n; i++){
2860 if (i==k)
2861 y[i] = 1.0;
2862 else
2863 y[i] = 0.0;
2864 for (j=0; j<=i-1; j++){
2865 y[i] = y[i] - a[i][j]*y[j];
2866 }
2867 }
2868
2869 /****************************************************
2870 Ux = y
2871 ****************************************************/
2872
2873 for (i=n; 0<=i; i--){
2874 x[i] = y[i];
2875 for (j=n; (i+1)<=j; j--){
2876 x[i] = x[i] - a[i][j]*x[j];
2877 }
2878 x[i] = x[i]/a[i][i];
2879 ia[i][k] = x[i];
2880 }
2881 }
2882
2883 for (i=0; i<=n; i++){
2884 for (j=0; j<=n; j++){
2885 a[i][j] = da[i][j];
2886 }
2887 }
2888 }
2889
2890 /***************************************************
2891 freeing of arrays:
2892
2893 x[List_YOUSO[39]]
2894 y[List_YOUSO[39]]
2895 da[List_YOUSO[39]][List_YOUSO[39]]
2896 ***************************************************/
2897
2898 free(x);
2899 free(y);
2900
2901 for (i=0; i<List_YOUSO[39]; i++){
2902 free(da[i]);
2903 }
2904 free(da);
2905 }
2906