1 /**********************************************************************
2 Hamiltonian_Cluster_NC.c:
3
4 Hamiltonian_Cluster_NC.c is a subroutine to make a Hamiltonian
5 matrix for cluster or molecular systems with non-collinear spin.
6
7 Log of Hamiltonian_Cluster_NC.c:
8
9 17/Dec/2003 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <math.h>
16 #include "openmx_common.h"
17 #include "mpi.h"
18
19
Hamiltonian_Cluster_NC(double ***** RH,double ***** IH,dcomplex ** H,int * MP)20 void Hamiltonian_Cluster_NC(double *****RH, double *****IH,
21 dcomplex **H, int *MP)
22 {
23 int i,j,k,q;
24 int MA_AN,GA_AN,LB_AN,GB_AN,AN;
25 int wanA,wanB,tnoA,tnoB,Anum,Bnum,NUM;
26 int num,tnum,num_orbitals,qmax;
27 int ID,myid,numprocs,tag=999;
28 int *My_NZeros;
29 int *is1,*ie1,*is2;
30 int *My_Matomnum,*order_GA;
31 double *H1,sum;
32
33 MPI_Status stat;
34 MPI_Request request;
35
36 /* MPI */
37
38 MPI_Comm_size(mpi_comm_level1,&numprocs);
39 MPI_Comm_rank(mpi_comm_level1,&myid);
40 MPI_Barrier(mpi_comm_level1);
41
42 /* allocation of arrays */
43
44 My_NZeros = (int*)malloc(sizeof(int)*numprocs);
45 My_Matomnum = (int*)malloc(sizeof(int)*numprocs);
46 is1 = (int*)malloc(sizeof(int)*numprocs);
47 ie1 = (int*)malloc(sizeof(int)*numprocs);
48 is2 = (int*)malloc(sizeof(int)*numprocs);
49 order_GA = (int*)malloc(sizeof(int)*(atomnum+2));
50
51 /* find my total number of non-zero elements in myid */
52
53 My_NZeros[myid] = 0;
54 for (MA_AN=1; MA_AN<=Matomnum; MA_AN++){
55 GA_AN = M2G[MA_AN];
56 wanA = WhatSpecies[GA_AN];
57 tnoA = Spe_Total_CNO[wanA];
58
59 num = 0;
60 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
61 GB_AN = natn[GA_AN][LB_AN];
62 wanB = WhatSpecies[GB_AN];
63 tnoB = Spe_Total_CNO[wanB];
64 num += tnoB;
65 }
66
67 My_NZeros[myid] += tnoA*num;
68 }
69
70 for (ID=0; ID<numprocs; ID++){
71 MPI_Bcast(&My_NZeros[ID],1,MPI_INT,ID,mpi_comm_level1);
72 }
73
74 tnum = 0;
75 for (ID=0; ID<numprocs; ID++){
76 tnum += My_NZeros[ID];
77 }
78
79 is1[0] = 0;
80 ie1[0] = My_NZeros[0] - 1;
81
82 for (ID=1; ID<numprocs; ID++){
83 is1[ID] = ie1[ID-1] + 1;
84 ie1[ID] = is1[ID] + My_NZeros[ID] - 1;
85 }
86
87 /* set is2 and order_GA */
88
89 My_Matomnum[myid] = Matomnum;
90 for (ID=0; ID<numprocs; ID++){
91 MPI_Bcast(&My_Matomnum[ID],1,MPI_INT,ID,mpi_comm_level1);
92 }
93
94 is2[0] = 1;
95 for (ID=1; ID<numprocs; ID++){
96 is2[ID] = is2[ID-1] + My_Matomnum[ID-1];
97 }
98
99 for (MA_AN=1; MA_AN<=Matomnum; MA_AN++){
100 order_GA[is2[myid]+MA_AN-1] = M2G[MA_AN];
101 }
102
103 for (ID=0; ID<numprocs; ID++){
104 MPI_Bcast(&order_GA[is2[ID]],My_Matomnum[ID],MPI_INT,ID,mpi_comm_level1);
105 }
106
107 /* set MP */
108
109 Anum = 1;
110 for (i=1; i<=atomnum; i++){
111 MP[i] = Anum;
112 wanA = WhatSpecies[i];
113 Anum = Anum + Spe_Total_CNO[wanA];
114 }
115 NUM = Anum - 1;
116 H[0][0].r = NUM;
117
118 /* set H1 */
119
120 H1 = (double*)malloc(sizeof(double)*(tnum+1));
121
122 /* initialize H */
123
124 for (i=1; i<=2*NUM; i++){
125 for (j=1; j<=2*NUM; j++){
126 H[i][j].r = 0.0;
127 H[i][j].i = 0.0;
128 }
129 }
130
131 /****************************************************************************
132 in case of SO_switch==0 && Hub_U_switch==0 && Constraint_NCS_switch==0
133 && Zeeman_NCS_switch==0 && Zeeman_NCO_switch==0
134
135 H[i ][j ].r = RH[0];
136 H[i ][j ].i = 0.0;
137 H[i+NUM][j+NUM].r = RH[1];
138 H[i+NUM][j+NUM].i = 0.0;
139 H[i ][j+NUM].r = RH[2];
140 H[i ][j+NUM].i = RH[3];
141
142 in case of SO_switch==1 or Hub_U_switch==1 or 1<=Constraint_NCS_switch
143 or Zeeman_NCS_switch==1 or Zeeman_NCO_switch==1
144
145 H[i ][j ].r = RH[0];
146 H[i ][j ].i = IH[0];
147 H[i+NUM][j+NUM].r = RH[1];
148 H[i+NUM][j+NUM].i = IH[1];
149 H[i ][j+NUM].r = RH[2];
150 H[i ][j+NUM].i = RH[3] + IH[2];
151 ****************************************************************************/
152
153 /* non-spin-orbit coupling and non-LDA+U */
154
155 if (SO_switch==0 && Hub_U_switch==0 && Constraint_NCS_switch==0
156 && Zeeman_NCS_switch==0 && Zeeman_NCO_switch==0)
157 qmax = 4;
158 else
159 qmax = 6;
160
161 for (q=0; q<qmax; q++){
162
163 k = is1[myid];
164 for (MA_AN=1; MA_AN<=Matomnum; MA_AN++){
165 GA_AN = M2G[MA_AN];
166 wanA = WhatSpecies[GA_AN];
167 tnoA = Spe_Total_CNO[wanA];
168 for (i=0; i<tnoA; i++){
169 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
170 GB_AN = natn[GA_AN][LB_AN];
171 wanB = WhatSpecies[GB_AN];
172 tnoB = Spe_Total_CNO[wanB];
173
174 /* non-spin-orbit coupling and non-LDA+U */
175
176 if (SO_switch==0 && Hub_U_switch==0 && Constraint_NCS_switch==0
177 && Zeeman_NCS_switch==0 && Zeeman_NCO_switch==0){
178 for (j=0; j<tnoB; j++){
179 H1[k] = RH[q][MA_AN][LB_AN][i][j];
180 k++;
181 }
182 }
183
184 /* spin-orbit coupling or LDA+U */
185
186 else{
187
188 switch(q){
189
190 case 0:
191
192 for (j=0; j<tnoB; j++){
193 H1[k] = RH[0][MA_AN][LB_AN][i][j];
194 k++;
195 }
196
197 break;
198
199 case 1:
200
201 for (j=0; j<tnoB; j++){
202 H1[k] = RH[1][MA_AN][LB_AN][i][j];
203 k++;
204 }
205
206 break;
207
208 case 2:
209
210 for (j=0; j<tnoB; j++){
211 H1[k] = RH[2][MA_AN][LB_AN][i][j];
212 k++;
213 }
214
215 break;
216
217 case 3:
218
219 for (j=0; j<tnoB; j++){
220 H1[k] = IH[0][MA_AN][LB_AN][i][j];
221 k++;
222 }
223
224 break;
225
226 case 4:
227
228 for (j=0; j<tnoB; j++){
229 H1[k] = IH[1][MA_AN][LB_AN][i][j];
230 k++;
231 }
232
233 break;
234
235 case 5:
236
237 for (j=0; j<tnoB; j++){
238 H1[k] = RH[3][MA_AN][LB_AN][i][j] + IH[2][MA_AN][LB_AN][i][j];
239 k++;
240 }
241
242 break;
243
244 } /* switch(q) */
245 } /* else */
246 }
247 }
248 }
249
250 /* MPI H1 */
251
252 for (ID=0; ID<numprocs; ID++){
253 k = is1[ID];
254 MPI_Bcast(&H1[k], My_NZeros[ID], MPI_DOUBLE, ID, mpi_comm_level1);
255 }
256
257 /* H1 -> H */
258
259 /* non-spin-orbit coupling and non-LDA+U */
260
261 if ( SO_switch==0 && Hub_U_switch==0 && Constraint_NCS_switch==0
262 && Zeeman_NCS_switch==0 && Zeeman_NCO_switch==0){
263
264 switch(q){
265
266 case 0:
267
268 k = 0;
269 for (AN=1; AN<=atomnum; AN++){
270 GA_AN = order_GA[AN];
271 wanA = WhatSpecies[GA_AN];
272 tnoA = Spe_Total_CNO[wanA];
273 Anum = MP[GA_AN];
274 for (i=0; i<tnoA; i++){
275 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
276 GB_AN = natn[GA_AN][LB_AN];
277 wanB = WhatSpecies[GB_AN];
278 tnoB = Spe_Total_CNO[wanB];
279 Bnum = MP[GB_AN];
280 for (j=0; j<tnoB; j++){
281 H[Anum+i][Bnum+j].r += H1[k];
282 k++;
283 }
284 }
285 }
286 }
287
288 break;
289
290 case 1:
291
292 k = 0;
293 for (AN=1; AN<=atomnum; AN++){
294 GA_AN = order_GA[AN];
295 wanA = WhatSpecies[GA_AN];
296 tnoA = Spe_Total_CNO[wanA];
297 Anum = MP[GA_AN];
298 for (i=0; i<tnoA; i++){
299 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
300 GB_AN = natn[GA_AN][LB_AN];
301 wanB = WhatSpecies[GB_AN];
302 tnoB = Spe_Total_CNO[wanB];
303 Bnum = MP[GB_AN];
304 for (j=0; j<tnoB; j++){
305 H[Anum+i+NUM][Bnum+j+NUM].r += H1[k];
306 k++;
307 }
308 }
309 }
310 }
311
312 break;
313
314 case 2:
315
316 k = 0;
317 for (AN=1; AN<=atomnum; AN++){
318 GA_AN = order_GA[AN];
319 wanA = WhatSpecies[GA_AN];
320 tnoA = Spe_Total_CNO[wanA];
321 Anum = MP[GA_AN];
322 for (i=0; i<tnoA; i++){
323 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
324 GB_AN = natn[GA_AN][LB_AN];
325 wanB = WhatSpecies[GB_AN];
326 tnoB = Spe_Total_CNO[wanB];
327 Bnum = MP[GB_AN];
328 for (j=0; j<tnoB; j++){
329 H[Anum+i][Bnum+j+NUM].r += H1[k];
330 H[Bnum+j+NUM][Anum+i].r += H1[k];
331 k++;
332 }
333 }
334 }
335 }
336
337 break;
338
339 case 3:
340
341 k = 0;
342 for (AN=1; AN<=atomnum; AN++){
343 GA_AN = order_GA[AN];
344 wanA = WhatSpecies[GA_AN];
345 tnoA = Spe_Total_CNO[wanA];
346 Anum = MP[GA_AN];
347 for (i=0; i<tnoA; i++){
348 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
349 GB_AN = natn[GA_AN][LB_AN];
350 wanB = WhatSpecies[GB_AN];
351 tnoB = Spe_Total_CNO[wanB];
352 Bnum = MP[GB_AN];
353 for (j=0; j<tnoB; j++){
354 H[Anum+i][Bnum+j+NUM].i += H1[k];
355 H[Bnum+j+NUM][Anum+i].i +=-H1[k];
356 k++;
357 }
358 }
359 }
360 }
361
362 break;
363
364 } /* switch(q) */
365 } /* if ( SO_switch==0 && Hub_U_switch==0 && Constraint_NCS_switch==0 && Zeeman_NCS_switch==0 && Zeeman_NCO_switch==0) */
366
367 /* in case of SO_switch==1 or Hub_U_switch==1 or 1<=Constraint_NCS_switch or Zeeman_NCS_switch==1 or Zeeman_NCO_switch==1 */
368
369 else{
370
371 switch(q){
372
373 case 0:
374
375 k = 0;
376 for (AN=1; AN<=atomnum; AN++){
377 GA_AN = order_GA[AN];
378 wanA = WhatSpecies[GA_AN];
379 tnoA = Spe_Total_CNO[wanA];
380 Anum = MP[GA_AN];
381 for (i=0; i<tnoA; i++){
382 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
383 GB_AN = natn[GA_AN][LB_AN];
384 wanB = WhatSpecies[GB_AN];
385 tnoB = Spe_Total_CNO[wanB];
386 Bnum = MP[GB_AN];
387 for (j=0; j<tnoB; j++){
388 H[Anum+i][Bnum+j].r += H1[k];
389 k++;
390 }
391 }
392 }
393 }
394
395 break;
396
397 case 1:
398
399 k = 0;
400 for (AN=1; AN<=atomnum; AN++){
401 GA_AN = order_GA[AN];
402 wanA = WhatSpecies[GA_AN];
403 tnoA = Spe_Total_CNO[wanA];
404 Anum = MP[GA_AN];
405 for (i=0; i<tnoA; i++){
406 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
407 GB_AN = natn[GA_AN][LB_AN];
408 wanB = WhatSpecies[GB_AN];
409 tnoB = Spe_Total_CNO[wanB];
410 Bnum = MP[GB_AN];
411 for (j=0; j<tnoB; j++){
412 H[Anum+i+NUM][Bnum+j+NUM].r += H1[k];
413 k++;
414 }
415 }
416 }
417 }
418
419 break;
420
421 case 2:
422
423 k = 0;
424 for (AN=1; AN<=atomnum; AN++){
425 GA_AN = order_GA[AN];
426 wanA = WhatSpecies[GA_AN];
427 tnoA = Spe_Total_CNO[wanA];
428 Anum = MP[GA_AN];
429 for (i=0; i<tnoA; i++){
430 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
431 GB_AN = natn[GA_AN][LB_AN];
432 wanB = WhatSpecies[GB_AN];
433 tnoB = Spe_Total_CNO[wanB];
434 Bnum = MP[GB_AN];
435 for (j=0; j<tnoB; j++){
436 H[Anum+i][Bnum+j+NUM].r += H1[k];
437 H[Bnum+j+NUM][Anum+i].r += H1[k];
438 k++;
439 }
440 }
441 }
442 }
443
444 break;
445
446 case 3:
447
448 k = 0;
449 for (AN=1; AN<=atomnum; AN++){
450 GA_AN = order_GA[AN];
451 wanA = WhatSpecies[GA_AN];
452 tnoA = Spe_Total_CNO[wanA];
453 Anum = MP[GA_AN];
454 for (i=0; i<tnoA; i++){
455 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
456 GB_AN = natn[GA_AN][LB_AN];
457 wanB = WhatSpecies[GB_AN];
458 tnoB = Spe_Total_CNO[wanB];
459 Bnum = MP[GB_AN];
460 for (j=0; j<tnoB; j++){
461 H[Anum+i][Bnum+j].i += H1[k];
462 k++;
463 }
464 }
465 }
466 }
467
468 break;
469
470 case 4:
471
472 k = 0;
473 for (AN=1; AN<=atomnum; AN++){
474 GA_AN = order_GA[AN];
475 wanA = WhatSpecies[GA_AN];
476 tnoA = Spe_Total_CNO[wanA];
477 Anum = MP[GA_AN];
478 for (i=0; i<tnoA; i++){
479 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
480 GB_AN = natn[GA_AN][LB_AN];
481 wanB = WhatSpecies[GB_AN];
482 tnoB = Spe_Total_CNO[wanB];
483 Bnum = MP[GB_AN];
484 for (j=0; j<tnoB; j++){
485 H[Anum+i+NUM][Bnum+j+NUM].i += H1[k];
486 k++;
487 }
488 }
489 }
490 }
491
492 break;
493
494 case 5:
495
496 k = 0;
497 for (AN=1; AN<=atomnum; AN++){
498 GA_AN = order_GA[AN];
499 wanA = WhatSpecies[GA_AN];
500 tnoA = Spe_Total_CNO[wanA];
501 Anum = MP[GA_AN];
502 for (i=0; i<tnoA; i++){
503 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
504 GB_AN = natn[GA_AN][LB_AN];
505 wanB = WhatSpecies[GB_AN];
506 tnoB = Spe_Total_CNO[wanB];
507 Bnum = MP[GB_AN];
508 for (j=0; j<tnoB; j++){
509 H[Anum+i][Bnum+j+NUM].i += H1[k];
510 H[Bnum+j+NUM][Anum+i].i +=-H1[k];
511 k++;
512 }
513 }
514 }
515 }
516
517 break;
518
519 } /* switch(q) */
520 } /* else */
521 } /* for (q=0;... */
522
523 /* freeing of arrays */
524
525 free(My_NZeros);
526 free(My_Matomnum);
527 free(is1);
528 free(ie1);
529 free(is2);
530 free(order_GA);
531 free(H1);
532 }
533
534
535