1 /**********************************************************************
2 Eff_Hub_Pot.c:
3
4 Eff_Hub_Pot.c is a subroutine to construct effective potential
5 for obtaining the effective Hubbard Hamiltonial in LDA+U calculation.
6
7 Log of Eff_Hub_Pot.c:
8
9 4/Aug/2004 Released by M.J.Han (--- MJ)
10 29/Nov/2004 Modified by T.Ozaki (AIST)
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <math.h>
16 #include <time.h>
17 #include "openmx_common.h"
18 #include "mpi.h"
19
20 static void H_U_onsite();
21 static void H_U_full(int SCF_iter, double ****OLP0);
22 static void H_U_dual(int SCF_iter, double ****OLP0);
23
Eff_Hub_Pot(int SCF_iter,double **** OLP0)24 void Eff_Hub_Pot(int SCF_iter, double ****OLP0)
25 {
26 /****************************************************
27 Important note:
28
29 In case of orbital optimization, the U potential
30 is applied to the primitive orbital.
31 ****************************************************/
32 int spin,a,b;
33 int l,fan,j,k,i,jg,kg,ig,wakg,waig,kl,il1,il2,m,n;
34 int Cwan,Hwan,Mul1,tno0,tno1,tno2;
35 int Mc_AN,Gc_AN,Mj_AN,Mi_AN,Mk_AN;
36 double sum;
37
38 /* in case of correction */
39
40 if ( Hub_U_switch==1
41 || 1<=Constraint_NCS_switch
42 || Zeeman_NCS_switch==1
43 || Zeeman_NCO_switch==1) {
44
45 /* on site */
46 if (Hub_U_occupation==0){
47 H_U_onsite();
48 }
49
50 /* full */
51 else if (Hub_U_occupation==1){
52 H_U_full(SCF_iter, OLP0);
53 }
54
55 /* dual */
56 else if (Hub_U_occupation==2){
57 H_U_dual(SCF_iter, OLP0);
58 }
59
60 } /* if */
61
62 /* in case of no correction */
63
64 else {
65
66 /* According to the report 44, iHNL0 is copied to iHNL. */
67
68 if (SpinP_switch==3 && SO_switch==1){
69
70 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
71
72 Gc_AN = M2G[Mc_AN];
73 Cwan = WhatSpecies[Gc_AN];
74 fan = FNAN[Gc_AN];
75 tno0 = Spe_Total_NO[Cwan];
76
77 for (j=0; j<=fan; j++){
78 jg = natn[Gc_AN][j];
79 Mj_AN = F_G2M[jg];
80 Hwan = WhatSpecies[jg];
81 tno1 = Spe_Total_NO[Hwan];
82
83 for (m=0; m<tno0; m++){
84 for (n=0; n<tno1; n++){
85
86 /* imaginary part */
87
88 iHNL[0][Mc_AN][j][m][n] = iHNL0[0][Mc_AN][j][m][n];
89 iHNL[1][Mc_AN][j][m][n] = iHNL0[1][Mc_AN][j][m][n];
90 iHNL[2][Mc_AN][j][m][n] = iHNL0[2][Mc_AN][j][m][n];
91 }
92 }
93 }
94 }
95 }
96
97 } /* else */
98
99 }
100
101
102
103
104
H_U_onsite()105 void H_U_onsite()
106 {
107 int spin,a,b;
108 int l,fan,j,k,i,jg,kg,ig,wakg,waig,kl,il1,il2,m,n;
109 int Cwan,Hwan,Mul1,tno0,tno1,tno2;
110 int Mc_AN,Gc_AN,Mj_AN,Mi_AN,Mk_AN;
111
112 /****************************************************
113 if (SpinP_switch!=3)
114
115 collinear case
116 ****************************************************/
117
118 if (SpinP_switch!=3){
119
120 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
121 Gc_AN = M2G[Mc_AN];
122 Cwan = WhatSpecies[Gc_AN];
123 fan = FNAN[Gc_AN];
124 tno0 = Spe_Total_NO[Cwan];
125
126 for (j=0; j<=fan; j++){
127 jg = natn[Gc_AN][j];
128 Mj_AN = F_G2M[jg];
129 Hwan = WhatSpecies[jg];
130 tno1 = Spe_Total_NO[Hwan];
131
132 if (j==0){
133 for (spin=0; spin<=SpinP_switch; spin++){
134 for (m=0; m<tno0; m++){
135 for (n=0; n<tno1; n++){
136 H_Hub[spin][Mc_AN][j][m][n] = v_eff[spin][Mc_AN][m][n];
137 }
138 }
139 }
140 }
141
142 else{
143 for (spin=0; spin<=SpinP_switch; spin++){
144 for (m=0; m<tno0; m++){
145 for (n=0; n<tno1; n++){
146 H_Hub[spin][Mc_AN][j][m][n] = 0.0;
147 }
148 }
149 }
150 }
151 }
152 }
153 }
154
155 /****************************************************
156 if (SpinP_switch==3)
157
158 spin non-collinear
159 ****************************************************/
160
161 else {
162
163 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
164 Gc_AN = M2G[Mc_AN];
165 Cwan = WhatSpecies[Gc_AN];
166 fan = FNAN[Gc_AN];
167 tno0 = Spe_Total_NO[Cwan];
168
169 for (j=0; j<=fan; j++){
170 jg = natn[Gc_AN][j];
171 Mj_AN = F_G2M[jg];
172 Hwan = WhatSpecies[jg];
173 tno1 = Spe_Total_NO[Hwan];
174
175 if (j==0){
176 for (m=0; m<tno0; m++){
177 for (n=0; n<tno1; n++){
178
179 /* real part */
180
181 H_Hub[0][Mc_AN][j][m][n] = NC_v_eff[0][0][Mc_AN][m][n].r;
182 H_Hub[1][Mc_AN][j][m][n] = NC_v_eff[1][1][Mc_AN][m][n].r;
183 H_Hub[2][Mc_AN][j][m][n] = NC_v_eff[0][1][Mc_AN][m][n].r;
184
185 /* imaginary part */
186
187 iHNL[0][Mc_AN][j][m][n] = iHNL0[0][Mc_AN][j][m][n] + F_U_flag*NC_v_eff[0][0][Mc_AN][m][n].i;
188 iHNL[1][Mc_AN][j][m][n] = iHNL0[1][Mc_AN][j][m][n] + F_U_flag*NC_v_eff[1][1][Mc_AN][m][n].i;
189 iHNL[2][Mc_AN][j][m][n] = iHNL0[2][Mc_AN][j][m][n] + F_U_flag*NC_v_eff[0][1][Mc_AN][m][n].i;
190 }
191 }
192 }
193
194 else{
195 for (m=0; m<tno0; m++){
196 for (n=0; n<tno1; n++){
197
198 /* real part */
199
200 H_Hub[0][Mc_AN][j][m][n] = 0.0;
201 H_Hub[1][Mc_AN][j][m][n] = 0.0;
202 H_Hub[2][Mc_AN][j][m][n] = 0.0;
203 }
204 }
205 }
206 }
207 }
208 }
209
210 }
211
212
H_U_full(int SCF_iter,double **** OLP0)213 void H_U_full(int SCF_iter, double ****OLP0)
214 {
215 int spin,a,b;
216 int l,fan,j,k,i,jg,kg,ig,wakg,waig,kl,il1,il2,m,n;
217 int Cwan,Hwan,Mul1,tno0,tno1,tno2;
218 int Mc_AN,Gc_AN,Mj_AN,Mi_AN,Mk_AN;
219 int num,h_AN,Gh_AN,size1,size2;
220 double *tmp_array;
221 double *tmp_array2;
222 int *Snd_S_Size,*Rcv_S_Size;
223 double sum,sum0;
224 double Resum00,Resum11,Resum01;
225 double Imsum00,Imsum11,Imsum01;
226 int numprocs,myid,ID,IDS,IDR,tag=999;
227
228 MPI_Status stat;
229 MPI_Request request;
230
231 /* MPI */
232 MPI_Comm_size(mpi_comm_level1,&numprocs);
233 MPI_Comm_rank(mpi_comm_level1,&myid);
234
235 /****************************************************
236 allocation of arrays:
237 ****************************************************/
238
239 Snd_S_Size = (int*)malloc(sizeof(int)*numprocs);
240 Rcv_S_Size = (int*)malloc(sizeof(int)*numprocs);
241
242 /****************************************************
243 MPI
244
245 OLP0
246 ****************************************************/
247
248 /***********************************
249 set data size
250 ************************************/
251
252 if (SCF_iter<=2){
253
254 for (ID=0; ID<numprocs; ID++){
255
256 IDS = (myid + ID) % numprocs;
257 IDR = (myid - ID + numprocs) % numprocs;
258
259 if (ID!=0){
260
261 tag = 999;
262
263 /* find data size to send block data */
264 if ((F_Snd_Num[IDS]+S_Snd_Num[IDS])!=0){
265 size1 = 0;
266 for (n=0; n<(F_Snd_Num[IDS]+S_Snd_Num[IDS]); n++){
267 Mc_AN = Snd_MAN[IDS][n];
268 Gc_AN = Snd_GAN[IDS][n];
269 Cwan = WhatSpecies[Gc_AN];
270 tno1 = Spe_Total_NO[Cwan];
271 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
272 Gh_AN = natn[Gc_AN][h_AN];
273 Hwan = WhatSpecies[Gh_AN];
274 tno2 = Spe_Total_NO[Hwan];
275 size1 += tno1*tno2;
276 }
277 }
278
279 Snd_S_Size[IDS] = size1;
280 MPI_Isend(&size1, 1, MPI_INT, IDS, tag, mpi_comm_level1, &request);
281 }
282 else{
283 Snd_S_Size[IDS] = 0;
284 }
285
286 /* receiving of size of data */
287
288 if ((F_Rcv_Num[IDR]+S_Rcv_Num[IDR])!=0){
289 MPI_Recv(&size2, 1, MPI_INT, IDR, tag, mpi_comm_level1, &stat);
290 Rcv_S_Size[IDR] = size2;
291 }
292 else{
293 Rcv_S_Size[IDR] = 0;
294 }
295
296 if ((F_Snd_Num[IDS]+S_Snd_Num[IDS])!=0) MPI_Wait(&request,&stat);
297 }
298 else{
299 Snd_S_Size[IDS] = 0;
300 Rcv_S_Size[IDR] = 0;
301 }
302 }
303
304 /***********************************
305 data transfer
306 ************************************/
307
308 tag = 999;
309 for (ID=0; ID<numprocs; ID++){
310
311 IDS = (myid + ID) % numprocs;
312 IDR = (myid - ID + numprocs) % numprocs;
313
314 if (ID!=0){
315
316 /*****************************
317 sending of data
318 *****************************/
319
320 if ((F_Snd_Num[IDS]+S_Snd_Num[IDS])!=0){
321
322 size1 = Snd_S_Size[IDS];
323
324 /* allocation of array */
325
326 tmp_array = (double*)malloc(sizeof(double)*size1);
327
328 /* multidimentional array to vector array */
329
330 num = 0;
331
332 for (n=0; n<(F_Snd_Num[IDS]+S_Snd_Num[IDS]); n++){
333 Mc_AN = Snd_MAN[IDS][n];
334 Gc_AN = Snd_GAN[IDS][n];
335 Cwan = WhatSpecies[Gc_AN];
336 tno1 = Spe_Total_NO[Cwan];
337 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
338 Gh_AN = natn[Gc_AN][h_AN];
339 Hwan = WhatSpecies[Gh_AN];
340 tno2 = Spe_Total_NO[Hwan];
341 for (i=0; i<tno1; i++){
342 for (j=0; j<tno2; j++){
343 tmp_array[num] = OLP0[Mc_AN][h_AN][i][j];
344 num++;
345 }
346 }
347 }
348 }
349
350 MPI_Isend(&tmp_array[0], size1, MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request);
351 }
352
353 /*****************************
354 receiving of block data
355 *****************************/
356
357 if ((F_Rcv_Num[IDR]+S_Rcv_Num[IDR])!=0){
358
359 size2 = Rcv_S_Size[IDR];
360
361 /* allocation of array */
362 tmp_array2 = (double*)malloc(sizeof(double)*size2);
363
364 MPI_Recv(&tmp_array2[0], size2, MPI_DOUBLE, IDR, tag, mpi_comm_level1, &stat);
365
366 num = 0;
367 Mc_AN = S_TopMAN[IDR] - 1; /* S_TopMAN should be used. */
368 for (n=0; n<(F_Rcv_Num[IDR]+S_Rcv_Num[IDR]); n++){
369 Mc_AN++;
370 Gc_AN = Rcv_GAN[IDR][n];
371 Cwan = WhatSpecies[Gc_AN];
372 tno1 = Spe_Total_NO[Cwan];
373
374 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
375 Gh_AN = natn[Gc_AN][h_AN];
376 Hwan = WhatSpecies[Gh_AN];
377 tno2 = Spe_Total_NO[Hwan];
378 for (i=0; i<tno1; i++){
379 for (j=0; j<tno2; j++){
380 OLP0[Mc_AN][h_AN][i][j] = tmp_array2[num];
381 num++;
382 }
383 }
384 }
385 }
386
387 /* freeing of array */
388 free(tmp_array2);
389 }
390
391 if ((F_Snd_Num[IDS]+S_Snd_Num[IDS])!=0){
392 MPI_Wait(&request,&stat);
393 free(tmp_array); /* freeing of array */
394 }
395 }
396 }
397 }
398
399 /****************************************************
400 make H_Hub
401 ****************************************************/
402
403 /****************************************************
404 if (SpinP_switch!=3)
405
406 collinear case
407 ****************************************************/
408
409 if (SpinP_switch!=3){
410
411 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
412 Gc_AN = M2G[Mc_AN];
413 Cwan = WhatSpecies[Gc_AN];
414 fan = FNAN[Gc_AN];
415 tno0 = Spe_Total_NO[Cwan];
416
417 for (j=0; j<=fan; j++){
418 jg = natn[Gc_AN][j];
419 Mj_AN = S_G2M[jg]; /* S_G2G should be used due to consistency with DC method */
420 Hwan = WhatSpecies[jg];
421 tno1 = Spe_Total_NO[Hwan];
422
423 for (spin=0; spin<=SpinP_switch; spin++){
424 for (m=0; m<tno0; m++){
425 for (n=0; n<tno1; n++){
426 H_Hub[spin][Mc_AN][j][m][n] = 0.0;
427 }
428 }
429 }
430
431 for (k=0; k<=fan; k++){
432 kg = natn[Gc_AN][k];
433 Mk_AN = F_G2M[kg]; /* F_G2G should be used */
434 wakg = WhatSpecies[kg];
435 kl = RMI1[Mc_AN][j][k];
436 tno2 = Spe_Total_NO[wakg];
437
438 if (0<=kl){
439
440 for (m=0; m<tno0; m++){
441 for (n=0; n<tno1; n++){
442 for (spin=0; spin<=SpinP_switch; spin++){
443
444 sum = 0.0;
445 for (a=0; a<tno2; a++){
446 for (b=0; b<tno2; b++){
447 sum += OLP0[Mc_AN][k][m][a] * v_eff[spin][Mk_AN][a][b] * OLP0[Mj_AN][kl][n][b];
448 }
449 }
450
451 if (3<=spin)
452 H_Hub[spin][Mc_AN][j][m][n] = 0.0;
453 else
454 H_Hub[spin][Mc_AN][j][m][n] += sum;
455 }
456 }
457 }
458 }
459 }
460 }
461 }
462 }
463
464 /****************************************************
465 if (SpinP_switch==3)
466
467 spin non-collinear
468 ****************************************************/
469
470 else{
471
472 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
473 Gc_AN = M2G[Mc_AN];
474 Cwan = WhatSpecies[Gc_AN];
475 fan = FNAN[Gc_AN];
476 tno0 = Spe_Total_NO[Cwan];
477
478 for (j=0; j<=fan; j++){
479 jg = natn[Gc_AN][j];
480 Mj_AN = S_G2M[jg]; /* S_G2G should be used due to consistency with DC method */
481 Hwan = WhatSpecies[jg];
482 tno1 = Spe_Total_NO[Hwan];
483
484 for (m=0; m<tno0; m++){
485 for (n=0; n<tno1; n++){
486
487 /* real part */
488
489 H_Hub[0][Mc_AN][j][m][n] = 0.0;
490 H_Hub[1][Mc_AN][j][m][n] = 0.0;
491 H_Hub[2][Mc_AN][j][m][n] = 0.0;
492
493 /* imaginary part */
494
495 iHNL[0][Mc_AN][j][m][n] = iHNL0[0][Mc_AN][j][m][n];
496 iHNL[1][Mc_AN][j][m][n] = iHNL0[1][Mc_AN][j][m][n];
497 iHNL[2][Mc_AN][j][m][n] = iHNL0[2][Mc_AN][j][m][n];
498 }
499 }
500
501 for (k=0; k<=fan; k++){
502 kg = natn[Gc_AN][k];
503 Mk_AN = F_G2M[kg]; /* F_G2G should be used */
504 wakg = WhatSpecies[kg];
505 kl = RMI1[Mc_AN][j][k];
506 tno2 = Spe_Total_NO[wakg];
507
508 if (0<=kl){
509
510 for (m=0; m<tno0; m++){
511 for (n=0; n<tno1; n++){
512
513 Resum00 = 0.0;
514 Resum11 = 0.0;
515 Resum01 = 0.0;
516
517 Imsum00 = 0.0;
518 Imsum11 = 0.0;
519 Imsum01 = 0.0;
520
521 for (a=0; a<tno2; a++){
522 for (b=0; b<tno2; b++){
523
524 Resum00 += OLP0[Mc_AN][k][m][a] * NC_v_eff[0][0][Mk_AN][a][b].r * OLP0[Mj_AN][kl][n][b];
525 Resum11 += OLP0[Mc_AN][k][m][a] * NC_v_eff[1][1][Mk_AN][a][b].r * OLP0[Mj_AN][kl][n][b];
526 Resum01 += OLP0[Mc_AN][k][m][a] * NC_v_eff[0][1][Mk_AN][a][b].r * OLP0[Mj_AN][kl][n][b];
527
528 Imsum00 += OLP0[Mc_AN][k][m][a] * NC_v_eff[0][0][Mk_AN][a][b].i * OLP0[Mj_AN][kl][n][b];
529 Imsum11 += OLP0[Mc_AN][k][m][a] * NC_v_eff[1][1][Mk_AN][a][b].i * OLP0[Mj_AN][kl][n][b];
530 Imsum01 += OLP0[Mc_AN][k][m][a] * NC_v_eff[0][1][Mk_AN][a][b].i * OLP0[Mj_AN][kl][n][b];
531 }
532 }
533
534 /* real part */
535
536 H_Hub[0][Mc_AN][j][m][n] += Resum00;
537 H_Hub[1][Mc_AN][j][m][n] += Resum11;
538 H_Hub[2][Mc_AN][j][m][n] += Resum01;
539
540 /* imaginary part */
541
542 iHNL[0][Mc_AN][j][m][n] += F_U_flag*Imsum00;
543 iHNL[1][Mc_AN][j][m][n] += F_U_flag*Imsum11;
544 iHNL[2][Mc_AN][j][m][n] += F_U_flag*Imsum01;
545
546 }
547 }
548 }
549 }
550 }
551 }
552 }
553
554 /****************************************************
555 freeing of arrays:
556 ****************************************************/
557
558 free(Snd_S_Size);
559 free(Rcv_S_Size);
560 }
561
562
563
564
565
566
H_U_dual(int SCF_iter,double **** OLP0)567 void H_U_dual(int SCF_iter, double ****OLP0)
568 {
569 int spin,a,b;
570 int l,fan,j,k,i,jg,kg,ig,wakg,waig,kl,il1,il2,m,n;
571 int Cwan,Hwan,Mul1,tno0,tno1,tno2;
572 int Mc_AN,Gc_AN,Mj_AN,Mi_AN,Mk_AN;
573 double tmp0;
574 double Resum00,Resum11,Resum01;
575 double Imsum00,Imsum11,Imsum01;
576
577 /****************************************************
578 if (SpinP_switch!=3)
579
580 collinear case
581 ****************************************************/
582
583 if (SpinP_switch!=3){
584
585 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
586
587 Gc_AN = M2G[Mc_AN];
588 Cwan = WhatSpecies[Gc_AN];
589 fan = FNAN[Gc_AN];
590 tno0 = Spe_Total_NO[Cwan];
591
592 for (j=0; j<=fan; j++){
593 jg = natn[Gc_AN][j];
594 Mj_AN = F_G2M[jg];
595 Hwan = WhatSpecies[jg];
596 tno1 = Spe_Total_NO[Hwan];
597
598 for (spin=0; spin<=SpinP_switch; spin++){
599 for (m=0; m<tno0; m++){
600 for (n=0; n<tno1; n++){
601
602 tmp0 = 0.0;
603
604 for (k=0; k<tno0; k++){
605 tmp0 += v_eff[spin][Mc_AN][m][k]*OLP0[Mc_AN][j][k][n];
606 }
607
608 for (k=0; k<tno1; k++){
609 tmp0 += v_eff[spin][Mj_AN][k][n]*OLP0[Mc_AN][j][m][k];
610 }
611
612 H_Hub[spin][Mc_AN][j][m][n] = 0.5*tmp0;
613 }
614 }
615 }
616 }
617 }
618 }
619
620 /****************************************************
621 if (SpinP_switch==3)
622
623 spin non-collinear
624 ****************************************************/
625
626 else{
627
628 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
629
630 Gc_AN = M2G[Mc_AN];
631 Cwan = WhatSpecies[Gc_AN];
632 fan = FNAN[Gc_AN];
633 tno0 = Spe_Total_NO[Cwan];
634
635 for (j=0; j<=fan; j++){
636 jg = natn[Gc_AN][j];
637 Mj_AN = F_G2M[jg];
638 Hwan = WhatSpecies[jg];
639 tno1 = Spe_Total_NO[Hwan];
640
641 for (m=0; m<tno0; m++){
642 for (n=0; n<tno1; n++){
643
644 Resum00 = 0.0;
645 Resum11 = 0.0;
646 Resum01 = 0.0;
647
648 Imsum00 = 0.0;
649 Imsum11 = 0.0;
650 Imsum01 = 0.0;
651
652 for (k=0; k<tno0; k++){
653
654 Resum00 += NC_v_eff[0][0][Mc_AN][m][k].r * OLP0[Mc_AN][j][k][n];
655 Resum11 += NC_v_eff[1][1][Mc_AN][m][k].r * OLP0[Mc_AN][j][k][n];
656 Resum01 += NC_v_eff[0][1][Mc_AN][m][k].r * OLP0[Mc_AN][j][k][n];
657
658 Imsum00 += NC_v_eff[0][0][Mc_AN][m][k].i * OLP0[Mc_AN][j][k][n];
659 Imsum11 += NC_v_eff[1][1][Mc_AN][m][k].i * OLP0[Mc_AN][j][k][n];
660 Imsum01 += NC_v_eff[0][1][Mc_AN][m][k].i * OLP0[Mc_AN][j][k][n];
661 }
662
663 for (k=0; k<tno1; k++){
664
665 Resum00 += NC_v_eff[0][0][Mj_AN][k][n].r * OLP0[Mc_AN][j][m][k];
666 Resum11 += NC_v_eff[1][1][Mj_AN][k][n].r * OLP0[Mc_AN][j][m][k];
667 Resum01 += NC_v_eff[0][1][Mj_AN][k][n].r * OLP0[Mc_AN][j][m][k];
668
669 Imsum00 += NC_v_eff[0][0][Mj_AN][k][n].i * OLP0[Mc_AN][j][m][k];
670 Imsum11 += NC_v_eff[1][1][Mj_AN][k][n].i * OLP0[Mc_AN][j][m][k];
671 Imsum01 += NC_v_eff[0][1][Mj_AN][k][n].i * OLP0[Mc_AN][j][m][k];
672 }
673
674 /* real part */
675
676 H_Hub[0][Mc_AN][j][m][n] = 0.5*Resum00;
677 H_Hub[1][Mc_AN][j][m][n] = 0.5*Resum11;
678 H_Hub[2][Mc_AN][j][m][n] = 0.5*Resum01;
679
680 /* imaginary part */
681
682 iHNL[0][Mc_AN][j][m][n] = iHNL0[0][Mc_AN][j][m][n] + 0.5*F_U_flag*Imsum00;
683 iHNL[1][Mc_AN][j][m][n] = iHNL0[1][Mc_AN][j][m][n] + 0.5*F_U_flag*Imsum11;
684 iHNL[2][Mc_AN][j][m][n] = iHNL0[2][Mc_AN][j][m][n] + 0.5*F_U_flag*Imsum01;
685 }
686 }
687 }
688 }
689 }
690
691 }
692
693
694