1 /**********************************************************************
2   Hamiltonian_Band_NC.c:
3 
4      Hamiltonian_Band_NC.c is a subroutine to make a spin non-collinear
5      Hamiltonian matrix for a periodic boundary system using Bloch theorem.
6 
7   Log of Hamiltonian_Band_NC.c:
8 
9      24/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 
Hamiltonian_Band_NC(int Host_ID1,double ***** RH,double ***** IH,dcomplex ** H,int * MP,double k1,double k2,double k3)19 void Hamiltonian_Band_NC(int Host_ID1,
20                          double *****RH, double *****IH,
21                          dcomplex **H, int *MP,
22                          double k1, double k2, double k3)
23 {
24   static int firsttime=1;
25   int i,j,k,wanA,wanB,tnoA,tnoB,Anum,Bnum;
26   int NUM,MA_AN,GA_AN,LB_AN,GB_AN;
27   int l1,l2,l3,Rn,n2;
28   double *tmp_array1,*tmp_array2;
29   double **H11r,**H11i;
30   double **H22r,**H22i;
31   double **H12r,**H12i;
32   double kRn,si,co,h;
33   int ID,myid,numprocs,tag=999;
34 
35   MPI_Status stat;
36   MPI_Request request;
37 
38   /* MPI */
39   MPI_Comm_size(mpi_comm_level1,&numprocs);
40   MPI_Comm_rank(mpi_comm_level1,&myid);
41   MPI_Barrier(mpi_comm_level1);
42 
43   /* set MP */
44   Anum = 1;
45   for (i=1; i<=atomnum; i++){
46     MP[i] = Anum;
47     wanA = WhatSpecies[i];
48     tnoA = Spe_Total_CNO[wanA];
49     Anum = Anum + tnoA;
50   }
51   NUM = Anum - 1;
52   n2 = NUM + 2;
53 
54   /*******************************************
55    allocation of H11r, H11i,
56                  H22r, H22i,
57                  H12r, H12i
58   *******************************************/
59 
60   H11r = (double**)malloc(sizeof(double*)*n2);
61   for (i=0; i<n2; i++){
62     H11r[i] = (double*)malloc(sizeof(double)*n2);
63   }
64 
65   H11i = (double**)malloc(sizeof(double*)*n2);
66   for (i=0; i<n2; i++){
67     H11i[i] = (double*)malloc(sizeof(double)*n2);
68   }
69 
70   H22r = (double**)malloc(sizeof(double*)*n2);
71   for (i=0; i<n2; i++){
72     H22r[i] = (double*)malloc(sizeof(double)*n2);
73   }
74 
75   H22i = (double**)malloc(sizeof(double*)*n2);
76   for (i=0; i<n2; i++){
77     H22i[i] = (double*)malloc(sizeof(double)*n2);
78   }
79 
80   H12r = (double**)malloc(sizeof(double*)*n2);
81   for (i=0; i<n2; i++){
82     H12r[i] = (double*)malloc(sizeof(double)*n2);
83   }
84 
85   H12i = (double**)malloc(sizeof(double*)*n2);
86   for (i=0; i<n2; i++){
87     H12i[i] = (double*)malloc(sizeof(double)*n2);
88   }
89 
90   /****************************************************
91    PrintMemory
92   ****************************************************/
93 
94   if (firsttime){
95     PrintMemory("Hamiltonian_Band: H11r",sizeof(double)*n2*n2,NULL);
96     PrintMemory("Hamiltonian_Band: H11i",sizeof(double)*n2*n2,NULL);
97     PrintMemory("Hamiltonian_Band: H22r",sizeof(double)*n2*n2,NULL);
98     PrintMemory("Hamiltonian_Band: H22i",sizeof(double)*n2*n2,NULL);
99     PrintMemory("Hamiltonian_Band: H12r",sizeof(double)*n2*n2,NULL);
100     PrintMemory("Hamiltonian_Band: H12i",sizeof(double)*n2*n2,NULL);
101   }
102 
103   /* for PrintMemory */
104   firsttime=0;
105 
106   /****************************************************
107                     set Hamiltonian
108   ****************************************************/
109 
110   H[0][0].r = 2.0*NUM;
111   for (i=1; i<=NUM; i++){
112     for (j=1; j<=NUM; j++){
113       H11r[i][j] = 0.0;
114       H11i[i][j] = 0.0;
115       H22r[i][j] = 0.0;
116       H22i[i][j] = 0.0;
117       H12r[i][j] = 0.0;
118       H12i[i][j] = 0.0;
119     }
120   }
121 
122   for (MA_AN=1; MA_AN<=Matomnum; MA_AN++){
123     GA_AN = M2G[MA_AN];
124     wanA = WhatSpecies[GA_AN];
125     tnoA = Spe_Total_CNO[wanA];
126     Anum = MP[GA_AN];
127 
128     for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
129       GB_AN = natn[GA_AN][LB_AN];
130       Rn = ncn[GA_AN][LB_AN];
131       wanB = WhatSpecies[GB_AN];
132       tnoB = Spe_Total_CNO[wanB];
133 
134       /*
135       kRn = k1*( (double)atv_ijk[Rn][1] + Cell_Gxyz[GB_AN][1] - Cell_Gxyz[GA_AN][1] )
136           + k2*( (double)atv_ijk[Rn][2] + Cell_Gxyz[GB_AN][2] - Cell_Gxyz[GA_AN][2] )
137           + k3*( (double)atv_ijk[Rn][3] + Cell_Gxyz[GB_AN][3] - Cell_Gxyz[GA_AN][3] );
138       */
139 
140       l1 = atv_ijk[Rn][1];
141       l2 = atv_ijk[Rn][2];
142       l3 = atv_ijk[Rn][3];
143       kRn = k1*(double)l1 + k2*(double)l2 + k3*(double)l3;
144 
145       si = sin(2.0*PI*kRn);
146       co = cos(2.0*PI*kRn);
147       Bnum = MP[GB_AN];
148 
149       /* non-spin-orbit coupling and non-LDA+U */
150       if (SO_switch==0 && Hub_U_switch==0 && Constraint_NCS_switch==0
151           && Zeeman_NCS_switch==0 && Zeeman_NCO_switch==0){
152 
153         for (i=0; i<=(tnoA-1); i++){
154 	  for (j=0; j<=(tnoB-1); j++){
155 	    H11r[Anum+i][Bnum+j] += co*RH[0][MA_AN][LB_AN][i][j];
156 	    H11i[Anum+i][Bnum+j] += si*RH[0][MA_AN][LB_AN][i][j];
157 	    H22r[Anum+i][Bnum+j] += co*RH[1][MA_AN][LB_AN][i][j];
158 	    H22i[Anum+i][Bnum+j] += si*RH[1][MA_AN][LB_AN][i][j];
159 	    H12r[Anum+i][Bnum+j] += co*RH[2][MA_AN][LB_AN][i][j] - si*RH[3][MA_AN][LB_AN][i][j];
160 	    H12i[Anum+i][Bnum+j] += si*RH[2][MA_AN][LB_AN][i][j] + co*RH[3][MA_AN][LB_AN][i][j];
161 	  }
162         }
163       }
164 
165       /* spin-orbit coupling or LDA+U */
166       else {
167         for (i=0; i<=(tnoA-1); i++){
168 	  for (j=0; j<=(tnoB-1); j++){
169 	    H11r[Anum+i][Bnum+j] += co*RH[0][MA_AN][LB_AN][i][j] - si*IH[0][MA_AN][LB_AN][i][j];
170 	    H11i[Anum+i][Bnum+j] += si*RH[0][MA_AN][LB_AN][i][j] + co*IH[0][MA_AN][LB_AN][i][j];
171 	    H22r[Anum+i][Bnum+j] += co*RH[1][MA_AN][LB_AN][i][j] - si*IH[1][MA_AN][LB_AN][i][j];
172 	    H22i[Anum+i][Bnum+j] += si*RH[1][MA_AN][LB_AN][i][j] + co*IH[1][MA_AN][LB_AN][i][j];
173 	    H12r[Anum+i][Bnum+j] += co*RH[2][MA_AN][LB_AN][i][j] - si*(RH[3][MA_AN][LB_AN][i][j]
174                                                                      + IH[2][MA_AN][LB_AN][i][j]);
175 	    H12i[Anum+i][Bnum+j] += si*RH[2][MA_AN][LB_AN][i][j] + co*(RH[3][MA_AN][LB_AN][i][j]
176                                                                      + IH[2][MA_AN][LB_AN][i][j]);
177 	  }
178         }
179       }
180 
181     }
182   }
183 
184   /******************************************************
185     MPI: H11r, H11i
186          H22r, H22i
187          H12r, H12i
188   ******************************************************/
189 
190   tmp_array1 = (double*)malloc(sizeof(double)*6*n2*List_YOUSO[7]);
191   tmp_array2 = (double*)malloc(sizeof(double)*6*n2*List_YOUSO[7]);
192 
193   if (myid!=Host_ID1){
194 
195     for (MA_AN=1; MA_AN<=Matomnum; MA_AN++){
196       GA_AN = M2G[MA_AN];
197       wanA = WhatSpecies[GA_AN];
198       tnoA = Spe_Total_CNO[wanA];
199       Anum = MP[GA_AN];
200 
201       k = 0;
202 
203       for (i=0; i<tnoA; i++){
204         for (j=0; j<n2; j++){
205           tmp_array1[k] = H11r[Anum+i][j];
206           k++;
207 	}
208       }
209 
210       for (i=0; i<tnoA; i++){
211         for (j=0; j<n2; j++){
212           tmp_array1[k] = H11i[Anum+i][j];
213           k++;
214 	}
215       }
216 
217       for (i=0; i<tnoA; i++){
218         for (j=0; j<n2; j++){
219           tmp_array1[k] = H22r[Anum+i][j];
220           k++;
221 	}
222       }
223 
224       for (i=0; i<tnoA; i++){
225         for (j=0; j<n2; j++){
226           tmp_array1[k] = H22i[Anum+i][j];
227           k++;
228 	}
229       }
230 
231       for (i=0; i<tnoA; i++){
232         for (j=0; j<n2; j++){
233           tmp_array1[k] = H12r[Anum+i][j];
234           k++;
235 	}
236       }
237 
238       for (i=0; i<tnoA; i++){
239         for (j=0; j<n2; j++){
240           tmp_array1[k] = H12i[Anum+i][j];
241           k++;
242 	}
243       }
244 
245       tag = 999;
246       MPI_Isend(&tmp_array1[0], 6*tnoA*n2, MPI_DOUBLE, Host_ID1,
247                 tag, mpi_comm_level1, &request);
248       MPI_Wait(&request,&stat);
249 
250     }
251 
252   }
253 
254   else{
255 
256     for (GA_AN=1; GA_AN<=atomnum; GA_AN++){
257       wanA = WhatSpecies[GA_AN];
258       tnoA = Spe_Total_CNO[wanA];
259       Anum = MP[GA_AN];
260       ID = G2ID[GA_AN];
261       if (ID!=Host_ID1){
262 
263         tag = 999;
264         MPI_Recv(&tmp_array2[0], 6*tnoA*n2, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
265 
266 	k = 0;
267 
268         for (i=0; i<tnoA; i++){
269           for (j=0; j<n2; j++){
270             H11r[Anum+i][j] = tmp_array2[k];
271 	    k++;
272 	  }
273 	}
274 
275         for (i=0; i<tnoA; i++){
276           for (j=0; j<n2; j++){
277             H11i[Anum+i][j] = tmp_array2[k];
278 	    k++;
279 	  }
280 	}
281 
282         for (i=0; i<tnoA; i++){
283           for (j=0; j<n2; j++){
284             H22r[Anum+i][j] = tmp_array2[k];
285 	    k++;
286 	  }
287 	}
288 
289         for (i=0; i<tnoA; i++){
290           for (j=0; j<n2; j++){
291             H22i[Anum+i][j] = tmp_array2[k];
292 	    k++;
293 	  }
294 	}
295 
296         for (i=0; i<tnoA; i++){
297           for (j=0; j<n2; j++){
298             H12r[Anum+i][j] = tmp_array2[k];
299 	    k++;
300 	  }
301 	}
302 
303         for (i=0; i<tnoA; i++){
304           for (j=0; j<n2; j++){
305             H12i[Anum+i][j] = tmp_array2[k];
306 	    k++;
307 	  }
308 	}
309 
310       }
311     }
312 
313   }
314 
315   free(tmp_array1);
316   free(tmp_array2);
317 
318   /******************************************************
319     the full complex matrix of H
320   ******************************************************/
321 
322   for (i=1; i<=NUM; i++){
323     for (j=1; j<=NUM; j++){
324       H[i    ][j    ].r =  H11r[i][j];
325       H[i    ][j    ].i =  H11i[i][j];
326       H[i+NUM][j+NUM].r =  H22r[i][j];
327       H[i+NUM][j+NUM].i =  H22i[i][j];
328       H[i    ][j+NUM].r =  H12r[i][j];
329       H[i    ][j+NUM].i =  H12i[i][j];
330       H[j+NUM][i    ].r =  H[i][j+NUM].r;
331       H[j+NUM][i    ].i = -H[i][j+NUM].i;
332     }
333   }
334 
335   /****************************************************
336                        free arrays
337   ****************************************************/
338 
339   for (i=0; i<n2; i++){
340     free(H11r[i]);
341   }
342   free(H11r);
343 
344   for (i=0; i<n2; i++){
345     free(H11i[i]);
346   }
347   free(H11i);
348 
349   for (i=0; i<n2; i++){
350     free(H22r[i]);
351   }
352   free(H22r);
353 
354   for (i=0; i<n2; i++){
355     free(H22i[i]);
356   }
357   free(H22i);
358 
359   for (i=0; i<n2; i++){
360     free(H12r[i]);
361   }
362   free(H12r);
363 
364   for (i=0; i<n2; i++){
365     free(H12i[i]);
366   }
367   free(H12i);
368 
369 }
370 
371 
372 
373 
374 
375