1 /**********************************************************************
2 EulerAngle_Spin.c:
3
4 EulerAngle_Spin.c is a subroutine to find the Euler angle of spin
5 orientation from the density matrix.
6
7 Log of EulerAngle_Spin.c:
8
9 15/Feb/2006 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <math.h>
15 #include "openmx_common.h"
16
17 #ifdef c_complex
18 #include <complex.h>
19 #endif
20
21 #define MIND 1.0e-8
22
23
EulerAngle_Spin(int quickcalc_flag,double Re11,double Re22,double Re12,double Im12,double Re21,double Im21,double Nup[2],double Ndown[2],double t[2],double p[2])24 void EulerAngle_Spin( int quickcalc_flag,
25 double Re11, double Re22,
26 double Re12, double Im12,
27 double Re21, double Im21,
28 double Nup[2], double Ndown[2],
29 double t[2], double p[2] )
30 {
31 double phi,theta,d1,d2,d3,d4;
32 double cop,sip,sit,cot,tmp,tmp1,tmp2,prod;
33 double mx,my,mz,tn,absm;
34 double S_coordinate[3];
35 dcomplex bunbo,bunsi;
36 dcomplex cd1,cd2,cd3,cd4,cNup,cNdown;
37
38 #ifdef c_complex
39 double complex ctmp1,ctheta,csit,ccot;
40 #endif
41
42
43 #ifdef c_complex
44
45 if (fabs(Re12)<1.0e-14){
46 phi = PI*90.0/180.0;
47 }
48 else{
49
50 bunbo.r = Re12 + Re21;
51 bunbo.i = Im12 + Im21;
52 bunsi.r = Re12 - Re21;
53 bunsi.i = Im12 - Im21;
54
55 tmp = -(bunsi.i*bunbo.r - bunsi.r*bunbo.i)/(bunbo.r*bunbo.r+bunbo.i*bunbo.i);
56 phi = atan(tmp);
57 }
58
59 cop = cos(phi);
60 sip = sin(phi);
61
62 if (fabs(Re11 - Re22)<1.0e-14){
63 ctheta = PI*90.0/180.0 + 0.0*I;
64 }
65 else {
66 tmp1 = (Re12*cop - Im12*sip + Re21*cop + Im21*sip)/(Re11 - Re22);
67 tmp2 = (Re12*sip + Im12*cop - Re21*sip + Im21*cop)/(Re11 - Re22);
68
69 ctmp1 = tmp1 + tmp2*I;
70 ctheta = catan(ctmp1);
71 }
72
73 csit = csin(ctheta);
74 ccot = ccos(ctheta);
75
76 cd1.r = 0.5*(Re11 + Re22);
77 cd1.i = 0.0;
78
79 cd2.r = 0.5*creal(ccot)*(Re11 - Re22);
80 cd2.i = 0.5*cimag(ccot)*(Re11 - Re22);
81
82 cd3.r = 0.5*( (Re12*creal(csit)-Im12*cimag(csit))*cop
83 -(Re12*cimag(csit)+Im12*creal(csit))*sip );
84
85 cd3.i = 0.5*( (Re12*creal(csit)-Im12*cimag(csit))*sip
86 +(Re12*cimag(csit)+Im12*creal(csit))*cop );
87
88 cd4.r = 0.5*( (Re21*creal(csit)-Im21*cimag(csit))*cop
89 +(Re21*cimag(csit)+Im21*creal(csit))*sip );
90
91 cd4.i = 0.5*(-(Re21*creal(csit)-Im21*cimag(csit))*sip
92 +(Re21*cimag(csit)+Im21*creal(csit))*cop );
93
94 cNup.r = cd1.r + cd2.r + cd3.r + cd4.r;
95 cNup.i = cd1.i + cd2.i + cd3.i + cd4.i;
96 cNdown.r = cd1.r - cd2.r - cd3.r - cd4.r;
97 cNdown.i = cd1.i - cd2.i - cd3.i - cd4.i;
98
99 Nup[0] = cNup.r;
100 Nup[1] = cNup.i;
101
102 Ndown[0] = cNdown.r;
103 Ndown[1] = cNdown.i;
104
105 t[0] = creal(ctheta);
106 t[1] = cimag(ctheta);
107
108 p[0] = phi;
109 p[1] = 0.0;
110
111 #else
112
113
114 /*
115 if (fabs(Re12)<MIND){
116
117 mx = 2.0*Re12;
118 my = -2.0*Im12;
119 phi = acos( mx/sqrt(mx*mx + my*my + MIND*MIND) );
120 }
121 else{
122 phi = atan(-Im12/Re12);
123 }
124
125 cop = cos(phi);
126 sip = sin(phi);
127
128 if (fabs(Re11 - Re22)<MIND){
129
130 mx = 2.0*Re12;
131 my = -2.0*Im12;
132 mz = Re11 - Re22;
133
134 theta = acos( mz/sqrt(mx*mx + my*my + mz*mz + MIND*MIND) );
135 }
136 else {
137 theta = atan( 2.0*(Re12*cop - Im12*sip) / (Re11 - Re22) );
138 }
139
140 sit = sin(theta);
141 cot = cos(theta);
142
143 d1 = 0.5*(Re11 + Re22);
144 d2 = 0.5*cot*(Re11 - Re22);
145 d3 = (Re12*cop - Im12*sip)*sit;
146
147 Nup[0] = d1 + d2 + d3;
148 Nup[1] = 0.0;
149
150 Ndown[0] = d1 - d2 - d3;
151 Ndown[1] = 0.0;
152
153 t[0] = theta;
154 t[1] = 0.0;
155
156 p[0] = phi;
157 p[1] = 0.0;
158 */
159
160
161
162 /*
163 printf("theta=%10.5f phi=%10.5f Re11=%10.5f Re22=%10.5f Re12=%10.5f Im12=%10.5f d1=%10.5f d2=%10.5f d3=%10.5f\n",
164 180.0*theta/PI,180.0*phi/PI,Re11,Re22,Re12,Im12,d1,d2,d3);
165 */
166
167
168
169 mx = 2.0*Re12;
170 my = -2.0*Im12;
171 mz = Re11 - Re22;
172
173 xyz2spherical(mx,my,mz,0.0,0.0,0.0,S_coordinate);
174
175 absm = S_coordinate[0];
176 theta = S_coordinate[1];
177 phi = S_coordinate[2];
178 tn = Re11 + Re22;
179
180 Nup[0] = 0.5*(tn + absm);
181 Nup[1] = 0.0;
182
183 Ndown[0] = 0.5*(tn - absm);
184 Ndown[1] = 0.0;
185
186 t[0] = theta;
187 t[1] = 0.0;
188
189 p[0] = phi;
190 p[1] = 0.0;
191
192 #endif
193
194 }
195