1 /**********************************************************************
2   Hamiltonian_Cluster.c:
3 
4      Hamiltonian_Cluster.c is a subroutine to make a Hamiltonian matrix
5      for cluster or molecular systems.
6 
7   Log of Hamiltonian_Cluster.c:
8 
9      22/Nov/2001  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 #ifndef scalapack
20 
Hamiltonian_Cluster_Hs(double **** RH,double * Hs,int * MP,int spin,int myworld1)21 void Hamiltonian_Cluster_Hs(double ****RH, double *Hs, int *MP, int spin, int myworld1)
22 {
23 }
24 
25 #endif
26 
27 #ifdef scalapack
28 
Hamiltonian_Cluster_Hs(double **** RH,double * Hs,int * MP,int spin,int myworld1)29 void Hamiltonian_Cluster_Hs(double ****RH, double *Hs, int *MP, int spin, int myworld1)
30 {
31   int i,j,k;
32   int MA_AN,GA_AN,LB_AN,GB_AN,AN;
33   int wanA,wanB,tnoA,tnoB,Anum,Bnum,NUM;
34   int num,tnum,num_orbitals;
35   int ID,myid,numprocs,tag=999;
36   int *My_NZeros;
37   int *is1,*ie1,*is2;
38   int *My_Matomnum,*order_GA;
39   double *H1,sum;
40 
41   MPI_Status stat;
42   MPI_Request request;
43 
44 
45   int ig,jg,il,jl,prow,pcol,brow,bcol;
46 
47   /* MPI */
48 
49   MPI_Comm_size(mpi_comm_level1,&numprocs);
50   MPI_Comm_rank(mpi_comm_level1,&myid);
51   MPI_Barrier(mpi_comm_level1);
52 
53   /* allocation of arrays */
54 
55   My_NZeros = (int*)malloc(sizeof(int)*numprocs);
56   My_Matomnum = (int*)malloc(sizeof(int)*numprocs);
57   is1 = (int*)malloc(sizeof(int)*numprocs);
58   ie1 = (int*)malloc(sizeof(int)*numprocs);
59   is2 = (int*)malloc(sizeof(int)*numprocs);
60   order_GA = (int*)malloc(sizeof(int)*(atomnum+2));
61 
62   /* find my total number of non-zero elements in myid */
63 
64   My_NZeros[myid] = 0;
65   for (MA_AN=1; MA_AN<=Matomnum; MA_AN++){
66     GA_AN = M2G[MA_AN];
67     wanA = WhatSpecies[GA_AN];
68     tnoA = Spe_Total_CNO[wanA];
69 
70     num = 0;
71     for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
72       GB_AN = natn[GA_AN][LB_AN];
73       wanB = WhatSpecies[GB_AN];
74       tnoB = Spe_Total_CNO[wanB];
75       num += tnoB;
76     }
77 
78     My_NZeros[myid] += tnoA*num;
79   }
80 
81   for (ID=0; ID<numprocs; ID++){
82     MPI_Bcast(&My_NZeros[ID],1,MPI_INT,ID,mpi_comm_level1);
83   }
84 
85   tnum = 0;
86   for (ID=0; ID<numprocs; ID++){
87     tnum += My_NZeros[ID];
88   }
89 
90   is1[0] = 0;
91   ie1[0] = My_NZeros[0] - 1;
92 
93   for (ID=1; ID<numprocs; ID++){
94     is1[ID] = ie1[ID-1] + 1;
95     ie1[ID] = is1[ID] + My_NZeros[ID] - 1;
96   }
97 
98   /* set is2 and order_GA */
99 
100   My_Matomnum[myid] = Matomnum;
101   for (ID=0; ID<numprocs; ID++){
102     MPI_Bcast(&My_Matomnum[ID],1,MPI_INT,ID,mpi_comm_level1);
103   }
104 
105   is2[0] = 1;
106   for (ID=1; ID<numprocs; ID++){
107     is2[ID] = is2[ID-1] + My_Matomnum[ID-1];
108   }
109 
110   for (MA_AN=1; MA_AN<=Matomnum; MA_AN++){
111     order_GA[is2[myid]+MA_AN-1] = M2G[MA_AN];
112   }
113 
114   for (ID=0; ID<numprocs; ID++){
115     MPI_Bcast(&order_GA[is2[ID]],My_Matomnum[ID],MPI_INT,ID,mpi_comm_level1);
116   }
117 
118   /* set MP */
119 
120   Anum = 1;
121   for (i=1; i<=atomnum; i++){
122     MP[i] = Anum;
123     wanA = WhatSpecies[i];
124     Anum += Spe_Total_CNO[wanA];
125   }
126   NUM = Anum - 1;
127 
128   /* set H1 */
129 
130   H1 = (double*)malloc(sizeof(double)*(tnum+1));
131 
132   k = is1[myid];
133   for (MA_AN=1; MA_AN<=Matomnum; MA_AN++){
134     GA_AN = M2G[MA_AN];
135     wanA = WhatSpecies[GA_AN];
136     tnoA = Spe_Total_CNO[wanA];
137     for (i=0; i<tnoA; i++){
138       for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
139         GB_AN = natn[GA_AN][LB_AN];
140         wanB = WhatSpecies[GB_AN];
141         tnoB = Spe_Total_CNO[wanB];
142         for (j=0; j<tnoB; j++){
143           H1[k] = RH[MA_AN][LB_AN][i][j];
144           k++;
145 	}
146       }
147     }
148   }
149 
150   /* MPI H1 */
151 
152   for (ID=0; ID<numprocs; ID++){
153     k = is1[ID];
154     MPI_Bcast(&H1[k], My_NZeros[ID], MPI_DOUBLE, ID, mpi_comm_level1);
155   }
156 
157   /* H1 -> H */
158 
159 
160   if(spin==myworld1){
161 
162     for(i=0;i<na_rows*na_cols;i++){
163       Hs[i] = 0.0;
164     }
165 
166 
167     k = 0;
168     for (AN=1; AN<=atomnum; AN++){
169       GA_AN = order_GA[AN];
170       wanA = WhatSpecies[GA_AN];
171       tnoA = Spe_Total_CNO[wanA];
172       Anum = MP[GA_AN];
173 
174       for (i=0; i<tnoA; i++){
175 
176 	for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
177 	  GB_AN = natn[GA_AN][LB_AN];
178 	  wanB = WhatSpecies[GB_AN];
179 	  tnoB = Spe_Total_CNO[wanB];
180 	  Bnum = MP[GB_AN];
181 
182 	  for (j=0; j<tnoB; j++){
183 	    ig = Anum+i;
184 	    jg = Bnum+j;
185 
186 	    brow = (ig-1)/nblk;
187 	    bcol = (jg-1)/nblk;
188 
189 	    prow = brow%np_rows;
190 	    pcol = bcol%np_cols;
191 
192 	    if(my_prow==prow && my_pcol==pcol){
193 	      il = (brow/np_rows+1)*nblk+1;
194 	      jl = (bcol/np_cols+1)*nblk+1;
195 	      if(((my_prow+np_rows)%np_rows) >= (brow%np_rows)){
196 		if(my_prow==prow){
197 		  il = il+(ig-1)%nblk;
198 		}
199 		il = il-nblk;
200 	      }
201 	      if(((my_pcol+np_cols)%np_cols) >= (bcol%np_cols)){
202 		if(my_pcol==pcol){
203 		  jl = jl+(jg-1)%nblk;
204 		}
205 		jl = jl-nblk;
206 	      }
207 	      Hs[(jl-1)*na_rows+il-1] += H1[k];
208 	    }
209 
210 	    k++;
211 	  }
212 	}
213       }
214     }
215 
216 
217   }
218 
219   /* freeing of arrays */
220 
221   free(My_NZeros);
222   free(My_Matomnum);
223   free(is1);
224   free(ie1);
225   free(is2);
226   free(order_GA);
227   free(H1);
228 }
229 
230 #endif
231