44 * is
55 * ( \vec{j}_s \cdot \nabla ) \vec{m}
66 *
7- * where j_s is the current vector in 2D
7+ * where j_s is the current vector in 3D
88 * For this computation we use the neighbours matrix to make easier the
99 * derivative calculation at the boundaries.
1010 * ** For now this only works with the CUBOID mesh
1717 * the derivatives are simply:
1818 * f(x + 1, y, z) - f(x - 1, y, z) / 2 dx
1919 * f(x, y + 1, z) - f(x, y - 1, z) / 2 dy
20+ * f(x, y, z + 1) - f(x, y, z - 1) / 2 dz
2021 *
2122 * where x+-1 is the neighbour to the left or right.
2223 *
3132 * i-th site and nn_x2 for the NN to the right. These variables are simply
3233 * the i-th index (in the corresponding cases) when there is a single NN
3334 *
34- * Similar for y
35+ * Similar for y and z
3536 *
3637 * n is the number of spins or lattice sites in the system
3738 *
5051 * Neighbouring sites where there is no material, has index -1
5152 *
5253 */
53- void compute_stt_field_c (double * spin , double * field , double * jx , double * jy ,
54- double dx , double dy , int * ngbs , int n ) {
54+ void compute_stt_field_c (double * spin , double * field , double * jx , double * jy , double * jz ,
55+ double dx , double dy , double dz , int * ngbs , int n ) {
5556
56- int i , j ;
57- float factor_x , factor_y ;
58- int nn_i ;
59- int nn_x1 , nn_x2 , nn_y1 , nn_y2 ;
60-
61- for (i = 0 ; i < 3 * n ; i ++ ) {
57+ //#pragma omp parallel for
58+ for (int i = 0 ; i < 3 * n ; i ++ ) {
6259 field [i ] = 0 ;
6360 }
64-
61+
62+ #pragma omp parallel for
6563 /* Iterate through every lattice site */
66- for (i = 0 ; i < n ; i ++ ){
64+ for (int i = 0 ; i < n ; i ++ ){
6765
6866 /* Starting index for the NNs of the i-th site
6967 * i+0, i+1, i+2, i+3 ... --> -x, +x, -y, +y ...
7068 */
71- nn_i = 6 * i ;
69+ int nn_i = 6 * i ;
70+ double factor_x , factor_y , factor_z ;
71+ int nn_x1 , nn_x2 , nn_y1 , nn_y2 , nn_z1 , nn_z2 ;
7272
7373 /* Here we distinguish if there are 2 NNs, no NN in the
7474 * -x direction, or no NN in the +x direction,
@@ -100,7 +100,7 @@ void compute_stt_field_c(double *spin, double *field, double *jx, double *jy,
100100 * This calculation is: jx[i] * d m[i] / dx
101101 * */
102102 if (factor_x ){
103- for (j = 0 ; j < 3 ; j ++ ){
103+ for (int j = 0 ; j < 3 ; j ++ ){
104104 field [3 * i + j ] += jx [i ] * (spin [3 * nn_x2 + j ]
105105 - spin [3 * nn_x1 + j ]) / (factor_x * dx );
106106 }
@@ -124,11 +124,37 @@ void compute_stt_field_c(double *spin, double *field, double *jx, double *jy,
124124 }
125125
126126 if (factor_y ){
127- for (j = 0 ; j < 3 ; j ++ ){
127+ for (int j = 0 ; j < 3 ; j ++ ){
128128 field [3 * i + j ] += jy [i ] * (spin [3 * nn_y2 + j ]
129129 - spin [3 * nn_y1 + j ]) / (factor_y * dy );
130130 }
131131 }
132+
133+
134+ // We do the same along the z direction
135+ if (ngbs [nn_i + 4 ] != -1 && ngbs [nn_i + 5 ] != -1 ) {
136+ factor_z = 2 ;
137+ nn_z1 = ngbs [nn_i + 4 ];
138+ nn_z2 = ngbs [nn_i + 5 ];
139+ } else if (ngbs [nn_i + 4 ] == -1 && ngbs [nn_i + 5 ] != -1 ){
140+ factor_z = 1 ;
141+ nn_z1 = ngbs [nn_i + 4 ];
142+ nn_z2 = i ;
143+ } else if (ngbs [nn_i + 4 ] == -1 && ngbs [nn_i + 5 ] != -1 ){
144+ factor_z = 1 ;
145+ nn_z1 = i ;
146+ nn_z2 = ngbs [nn_i + 5 ];
147+ } else {
148+ factor_z = 0 ;
149+ }
150+
151+ if (factor_z ){
152+ for (int j = 0 ; j < 3 ; j ++ ){
153+ field [3 * i + j ] += jz [i ] * (spin [3 * nn_z2 + j ]
154+ - spin [3 * nn_z1 + j ]) / (factor_z * dz );
155+ }
156+ }
157+
132158 }
133159}
134160
0 commit comments