1 /**********************************************************************
2 TRAN_Calc_GridBound.c:
3
4 TRAN_Calc_GridBound.c is a subroutine to find the grid boundary region.
5
6 output: TRAN_region[] is changed
7 TRAN_grid_bound[2];
8
9 The boundary regions are defined by [0:TRAN_grid_bound[0]] and [TRAN_grid_bound[1]: Ngrid1-1]
10
11
12 Log of TRAN_Calc_GridBound.c:
13
14 24/July/2008 Released by H.Kino and T.Ozaki
15
16 ***********************************************************************/
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <string.h>
21 #include <math.h>
22 #include <mpi.h>
23 #include "tran_variables.h"
24 #include "tran_prototypes.h"
25
26
27 #ifdef MAX
28 #undef MAX
29 #endif
30 #define MAX(a,b) (((a)>(b))? (a):(b))
31
32 #ifdef MIN
33 #undef MIN
34 #endif
35 #define MIN(a,b) (((a)<(b))? (a):(b))
36
37
38 static void Cross_Product(double a[4], double b[4], double c[4]);
39 static double Dot_Product(double a[4], double b[4]);
40
41
TRAN_Calc_GridBound(MPI_Comm mpi_comm_level1,int atomnum,int * WhatSpecies,double * Spe_Atom_Cut1,int Ngrid1,double * Grid_Origin,double ** Gxyz,double tv[4][4],double gtv[4][4],double rgtv[4][4],double Left_tv[4][4],double Right_tv[4][4])42 void TRAN_Calc_GridBound(MPI_Comm mpi_comm_level1,
43 int atomnum,
44 int *WhatSpecies,
45 double *Spe_Atom_Cut1,
46 int Ngrid1,
47 double *Grid_Origin,
48 double **Gxyz,
49 double tv[4][4],
50 double gtv[4][4],
51 double rgtv[4][4],
52 double Left_tv[4][4],
53 double Right_tv[4][4])
54 {
55 int ct_AN,ct_AN2,wanA,wanB;
56 int i,j,po;
57 double r,rcutA,rcutB;
58 double dx,dy,dz;
59 double MinV,MaxV;
60 int MaxGridNum;
61 double R[4],Cxyz[4];
62 double b[4],c[4],v[4];
63 double coef,Vec0,Vec1;
64 int myid,numprocs;
65
66 /* MPI */
67 MPI_Comm_size(mpi_comm_level1,&numprocs);
68 MPI_Comm_rank(mpi_comm_level1,&myid);
69
70 for (i=1; i<=3; i++) { R[i] = tv[1][i] - Right_tv[1][i]; }
71
72 if (myid==Host_ID){
73 printf("\n\n<TRAN_Calc_GridBound>\n\n");
74 }
75
76 /******************************************************
77 find the unit vector perpendicular to the bc-plane
78 ******************************************************/
79
80 b[1] = tv[2][1];
81 b[2] = tv[2][2];
82 b[3] = tv[2][3];
83
84 c[1] = tv[3][1];
85 c[2] = tv[3][2];
86 c[3] = tv[3][3];
87
88 Cross_Product(b,c,v);
89 coef = 1.0/sqrt(fabs( Dot_Product(v,v) ));
90
91 v[1] = coef*v[1];
92 v[2] = coef*v[2];
93 v[3] = coef*v[3];
94
95 /**********************************************************
96 left side
97 ***********************************************************/
98
99 MinV = 1.0e+10;
100 MaxV = -1.0e+10;
101
102 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
103
104 if ( TRAN_region[ct_AN]==2 ) { /* left region */
105
106 wanA = WhatSpecies[ct_AN];
107 rcutA = Spe_Atom_Cut1[wanA];
108
109 /* find Vec0 and Vec1 */
110
111 Cxyz[1] = Gxyz[ct_AN][1] + rcutA*v[1] - Grid_Origin[1];
112 Cxyz[2] = Gxyz[ct_AN][2] + rcutA*v[2] - Grid_Origin[2];
113 Cxyz[3] = Gxyz[ct_AN][3] + rcutA*v[3] - Grid_Origin[3];
114
115 Vec0 = Dot_Product(Cxyz,rgtv[1])*0.5/PI;
116
117 Cxyz[1] = Gxyz[ct_AN][1] - rcutA*v[1] - Grid_Origin[1];
118 Cxyz[2] = Gxyz[ct_AN][2] - rcutA*v[2] - Grid_Origin[2];
119 Cxyz[3] = Gxyz[ct_AN][3] - rcutA*v[3] - Grid_Origin[3];
120
121 Vec1 = Dot_Product(Cxyz,rgtv[1])*0.5/PI;
122
123 if (Vec0<MinV) MinV = Vec0;
124 if (Vec1<MinV) MinV = Vec1;
125 if (MaxV<Vec0) MaxV = Vec0;
126 if (MaxV<Vec1) MaxV = Vec1;
127
128 /* set "12" atoms having non-zero overlap with atoms in the left lead */
129
130 ct_AN2 = 1;
131 po = 0;
132
133 do {
134
135 if ( TRAN_region[ct_AN2]==2 ) { /* left region */
136
137 wanB = WhatSpecies[ct_AN2];
138 rcutB = Spe_Atom_Cut1[wanB];
139
140 dx = Gxyz[ct_AN][1] - (Gxyz[ct_AN2][1] - Left_tv[1][1]);
141 dy = Gxyz[ct_AN][2] - (Gxyz[ct_AN2][2] - Left_tv[1][2]);
142 dz = Gxyz[ct_AN][3] - (Gxyz[ct_AN2][3] - Left_tv[1][3]);
143
144 r = sqrt(dx*dx+dy*dy+dz*dz);
145
146 if (r<(rcutA+rcutB) ){
147 TRAN_region[ct_AN] += 10;
148 po = 1;
149 }
150 }
151
152 ct_AN2++;
153
154 } while (po==0 && ct_AN2<=atomnum);
155 }
156 }
157
158 /***************************************************
159 TRAN_grid_bound[0] is the maximum grid point that
160 basis functions of atoms belonging to the left lead
161 can ovelap.
162 ***************************************************/
163
164 Vec0 = fabs(Dot_Product(Left_tv[1],rgtv[1])*0.5/PI);
165 TRAN_grid_bound[0] = (int)(MaxV - Vec0) + 1;
166
167 /*
168 if (myid==Host_ID){
169 printf("Left MinV=%15.12f MaxV=%15.12f\n",MinV,MaxV);
170 printf("grid_bound (left) =%d\n",TRAN_grid_bound[0]);
171 }
172 */
173
174 /**********************************************************
175 right side
176 ***********************************************************/
177
178 MinV = 1.0e+10;
179 MaxV = -1.0e+10;
180
181 MaxGridNum = Dot_Product(tv[1],rgtv[1])*0.5/PI;
182
183 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
184
185 if ( TRAN_region[ct_AN]==3 ) { /* right region */
186
187 wanA = WhatSpecies[ct_AN];
188 rcutA = Spe_Atom_Cut1[wanA];
189
190 /* find Vec0 and Vec1 */
191
192 Cxyz[1] = Gxyz[ct_AN][1] + rcutA*v[1] - Grid_Origin[1];
193 Cxyz[2] = Gxyz[ct_AN][2] + rcutA*v[2] - Grid_Origin[2];
194 Cxyz[3] = Gxyz[ct_AN][3] + rcutA*v[3] - Grid_Origin[3];
195
196 Vec0 = Dot_Product(Cxyz,rgtv[1])*0.5/PI;
197
198 Cxyz[1] = Gxyz[ct_AN][1] - rcutA*v[1] - Grid_Origin[1];
199 Cxyz[2] = Gxyz[ct_AN][2] - rcutA*v[2] - Grid_Origin[2];
200 Cxyz[3] = Gxyz[ct_AN][3] - rcutA*v[3] - Grid_Origin[3];
201
202 Vec1 = Dot_Product(Cxyz,rgtv[1])*0.5/PI;
203
204 if (Vec0<MinV) MinV = Vec0;
205 if (Vec1<MinV) MinV = Vec1;
206 if (MaxV<Vec0) MaxV = Vec0;
207 if (MaxV<Vec1) MaxV = Vec1;
208
209 /* set "13" atoms having non-zero overlap with atoms in the right lead */
210
211 ct_AN2 = 1;
212 po = 0;
213
214 do {
215
216 if ( TRAN_region[ct_AN2]==3 ) { /* right region */
217
218 wanB = WhatSpecies[ct_AN2];
219 rcutB = Spe_Atom_Cut1[wanB];
220
221 dx = Gxyz[ct_AN][1] - (Gxyz[ct_AN2][1] + Right_tv[1][1]);
222 dy = Gxyz[ct_AN][2] - (Gxyz[ct_AN2][2] + Right_tv[1][2]);
223 dz = Gxyz[ct_AN][3] - (Gxyz[ct_AN2][3] + Right_tv[1][3]);
224
225 r = sqrt(dx*dx+dy*dy+dz*dz);
226
227 if (r<(rcutA+rcutB) ){
228 TRAN_region[ct_AN] += 10;
229 po = 1;
230 }
231 }
232
233 ct_AN2++;
234
235 } while (po==0 && ct_AN2<=atomnum);
236
237 }
238 }
239
240 /***************************************************
241 TRAN_grid_bound[1] is the minimum grid point that
242 basis functions of atoms belonging to the right lead
243 can ovelap.
244 ***************************************************/
245
246 Vec0 = fabs(Dot_Product(Right_tv[1],rgtv[1])*0.5/PI);
247 TRAN_grid_bound[1]= (int)(MinV + Vec0);
248
249 /*
250 if (myid==Host_ID){
251 printf("Right MinV=%15.12f MaxV=%15.12f\n",MinV,MaxV);
252 printf("grid_bound (right) =%d\n",TRAN_grid_bound[1]);
253 }
254 */
255
256 if (myid==Host_ID){
257
258 printf("\n*******************************************************\n"); fflush(stdout);
259 printf("The extended cell consists of Left0-Center-Right0.\n"); fflush(stdout);
260 printf("The cells of left and right reads are connected as.\n"); fflush(stdout);
261 printf("...|Left2|Left1|Left0-Center-Right0|Right1|Right2...\n\n"); fflush(stdout);
262 printf("Each atom in the extended cell is assigned as follows:\n"); fflush(stdout);
263 printf("where '12' and '2' mean that they are in 'Left0', and\n"); fflush(stdout);
264 printf("'12' has overlap with atoms in the Left1,\n"); fflush(stdout);
265 printf("and '13' and '3' mean that they are in 'Right0', and\n"); fflush(stdout);
266 printf("'13' has overlap with atoms in the 'Right1', and also\n"); fflush(stdout);
267 printf("'1' means atom in the 'Center'.\n"); fflush(stdout);
268 printf("********************************************************\n\n"); fflush(stdout);
269
270 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
271 if (ct_AN<10)
272 printf("Atom%d =%3d ",ct_AN,TRAN_region[ct_AN]);
273 else if (ct_AN<100)
274 printf("Atom%d =%3d ",ct_AN,TRAN_region[ct_AN]);
275 else if (ct_AN<1000)
276 printf("Atom%d=%3d ",ct_AN,TRAN_region[ct_AN]);
277
278 if (ct_AN%7==0) printf("\n");
279 }
280 printf("\n\n");
281 }
282
283 }
284
285
286
287
Cross_Product(double a[4],double b[4],double c[4])288 void Cross_Product(double a[4], double b[4], double c[4])
289 {
290 c[1] = a[2]*b[3] - a[3]*b[2];
291 c[2] = a[3]*b[1] - a[1]*b[3];
292 c[3] = a[1]*b[2] - a[2]*b[1];
293 }
294
Dot_Product(double a[4],double b[4])295 double Dot_Product(double a[4], double b[4])
296 {
297 double sum;
298 sum = a[1]*b[1] + a[2]*b[2] + a[3]*b[3];
299 return sum;
300 }
301