1 /**********************************************************************
2 Simple_Mixing_DM.c:
3
4 Simple_Mixing_DM.c is a subroutine to achieve self-consistent
5 field using the simple mixing method.
6
7 Log of Simple_Mixing_DM.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
Simple_Mixing_DM(int Change_switch,double Mix_wgt,double ***** CDM,double ***** PDM,double ***** P2DM,double ***** iCDM,double ***** iPDM,double ***** iP2DM,double ***** RDM,double ***** iRDM,EXX_t * exx,dcomplex **** exx_CDM,dcomplex **** exx_PDM,dcomplex **** exx_P2DM)20 void Simple_Mixing_DM(int Change_switch,
21 double Mix_wgt,
22 double *****CDM,
23 double *****PDM,
24 double *****P2DM,
25 double *****iCDM,
26 double *****iPDM,
27 double *****iP2DM,
28 double *****RDM,
29 double *****iRDM,
30 /*---------- modified by TOYODA 18/JAN/2010 */
31 EXX_t *exx,
32 dcomplex ****exx_CDM,
33 dcomplex ****exx_PDM,
34 dcomplex ****exx_P2DM
35 /*---------- until here */
36 )
37 {
38 int ian,jan,Mc_AN,Gc_AN;
39 int h_AN,Gh_AN,m,n,spin,i;
40 double Mix_wgt1,Mix_wgt2,Norm,My_Norm,tmp0;
41 double Min_Weight,Max_Weight;
42 double nc_weight[4];
43 int numprocs,myid,ID;
44
45 /* MPI */
46 MPI_Comm_size(mpi_comm_level1,&numprocs);
47 MPI_Comm_rank(mpi_comm_level1,&myid);
48
49 /* start... */
50
51 Min_Weight = Min_Mixing_weight;
52 if (SCF_RENZOKU==-1){
53 Max_Weight = Max_Mixing_weight;
54 Max_Mixing_weight2 = Max_Mixing_weight;
55 }
56 else if (SCF_RENZOKU==1000){ /* past 3 */
57 Max_Mixing_weight2 = 2.0*Max_Mixing_weight2;
58 if (0.7<Max_Mixing_weight2) Max_Mixing_weight2 = 0.7;
59 Max_Weight = Max_Mixing_weight2;
60 SCF_RENZOKU = 0;
61 }
62 else{
63 Max_Weight = Max_Mixing_weight2;
64 }
65
66 /****************************************************
67 NormRD
68 ****************************************************/
69
70 My_Norm = 0.0;
71
72 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
73 Gc_AN = M2G[Mc_AN];
74 ian = Spe_Total_CNO[WhatSpecies[Gc_AN]];
75 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
76 Gh_AN = natn[Gc_AN][h_AN];
77 jan = Spe_Total_CNO[WhatSpecies[Gh_AN]];
78
79 for (spin=0; spin<=SpinP_switch; spin++){
80 for (m=0; m<ian; m++){
81 for (n=0; n<jan; n++){
82
83 RDM[spin][Mc_AN][h_AN][m][n] = CDM[spin][Mc_AN][h_AN][m][n]
84 -PDM[spin][Mc_AN][h_AN][m][n];
85 My_Norm += RDM[spin][Mc_AN][h_AN][m][n]*RDM[spin][Mc_AN][h_AN][m][n];
86 }
87 }
88 }
89
90 if (SpinP_switch==3 && ( SO_switch==1 || Hub_U_switch==1 || 1<=Constraint_NCS_switch
91 || Zeeman_NCS_switch==1 || Zeeman_NCO_switch==1 )){
92
93 for (spin=0; spin<2; spin++){
94 for (m=0; m<ian; m++){
95 for (n=0; n<jan; n++){
96
97 iRDM[spin][Mc_AN][h_AN][m][n] = iCDM[spin][Mc_AN][h_AN][m][n]
98 -iPDM[spin][Mc_AN][h_AN][m][n];
99 My_Norm += iRDM[spin][Mc_AN][h_AN][m][n]*iRDM[spin][Mc_AN][h_AN][m][n];
100 }
101 }
102 }
103 }
104 }
105 }
106
107 /****************************************************
108 MPI:
109
110 My_Norm
111 ****************************************************/
112
113 MPI_Allreduce(&My_Norm, &Norm, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
114
115 /****************************************************
116 find an optimum mixing weight
117 ****************************************************/
118
119 for (i=4; 1<=i; i--){
120 NormRD[i] = NormRD[i-1];
121 History_Uele[i] = History_Uele[i-1];
122 }
123 NormRD[0] = Norm;
124 History_Uele[0] = Uele;
125
126 if (Change_switch==1){
127
128 if ( (int)sgn(History_Uele[0]-History_Uele[1])
129 ==(int)sgn(History_Uele[1]-History_Uele[2])
130 && NormRD[0]<NormRD[1]){
131
132 /* tmp0 = 1.6*Mixing_weight; */
133
134 tmp0 = NormRD[1]/(largest(NormRD[1]-NormRD[0],10e-10))*Mixing_weight;
135
136 if (tmp0<Max_Weight){
137 if (Min_Weight<tmp0){
138 Mixing_weight = tmp0;
139 }
140 else{
141 Mixing_weight = Min_Weight;
142 }
143 }
144 else{
145 Mixing_weight = Max_Weight;
146 SCF_RENZOKU++;
147 }
148 }
149
150 else if ( (int)sgn(History_Uele[0]-History_Uele[1])
151 ==(int)sgn(History_Uele[1]-History_Uele[2])
152 && NormRD[1]<NormRD[0]){
153
154 tmp0 = NormRD[1]/(largest(NormRD[1]+NormRD[0],10e-10))*Mixing_weight;
155
156 /* tmp0 = Mixing_weight/1.6; */
157
158 if (tmp0<Max_Weight){
159 if (Min_Weight<tmp0)
160 Mixing_weight = tmp0;
161 else
162 Mixing_weight = Min_Weight;
163 }
164 else
165 Mixing_weight = Max_Weight;
166
167 SCF_RENZOKU = -1;
168 }
169
170 else if ( (int)sgn(History_Uele[0]-History_Uele[1])
171 !=(int)sgn(History_Uele[1]-History_Uele[2])
172 && NormRD[0]<NormRD[1]){
173
174 /* tmp0 = Mixing_weight*1.2; */
175
176 tmp0 = NormRD[1]/(largest(NormRD[1]-NormRD[0],10e-10))*Mixing_weight;
177
178 if (tmp0<Max_Weight){
179 if (Min_Weight<tmp0)
180 Mixing_weight = tmp0;
181 else
182 Mixing_weight = Min_Weight;
183 }
184 else{
185 Mixing_weight = Max_Weight;
186 SCF_RENZOKU++;
187 }
188 }
189
190 else if ( (int)sgn(History_Uele[0]-History_Uele[1])
191 !=(int)sgn(History_Uele[1]-History_Uele[2])
192 && NormRD[1]<NormRD[0]){
193
194 /* tmp0 = Mixing_weight/2.0; */
195
196 tmp0 = NormRD[1]/(largest(NormRD[1]+NormRD[0],10e-10))*Mixing_weight;
197
198 if (tmp0<Max_Weight){
199 if (Min_Weight<tmp0)
200 Mixing_weight = tmp0;
201 else
202 Mixing_weight = Min_Weight;
203 }
204 else
205 Mixing_weight = Max_Weight;
206
207 SCF_RENZOKU = -1;
208 }
209
210 Mix_wgt = Mixing_weight;
211 }
212
213 /****************************************************
214 Mixing
215 ****************************************************/
216
217 nc_weight[0] = 1.0;
218 nc_weight[1] = 1.0;
219 nc_weight[2] = 1.0;
220 nc_weight[3] = 1.0;
221
222 for (spin=0; spin<=SpinP_switch; spin++){
223
224 Mix_wgt1 = nc_weight[spin]*Mix_wgt;
225 Mix_wgt2 = 1.0 - Mix_wgt1;
226
227 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
228 Gc_AN = M2G[Mc_AN];
229 ian = Spe_Total_CNO[WhatSpecies[Gc_AN]];
230 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
231 Gh_AN = natn[Gc_AN][h_AN];
232 jan = Spe_Total_CNO[WhatSpecies[Gh_AN]];
233
234 for (m=0; m<ian; m++){
235 for (n=0; n<jan; n++){
236
237 CDM[spin][Mc_AN][h_AN][m][n] = Mix_wgt1 * CDM[spin][Mc_AN][h_AN][m][n]
238 + Mix_wgt2 * PDM[spin][Mc_AN][h_AN][m][n];
239
240 P2DM[spin][Mc_AN][h_AN][m][n] = PDM[spin][Mc_AN][h_AN][m][n];
241 PDM[spin][Mc_AN][h_AN][m][n] = CDM[spin][Mc_AN][h_AN][m][n];
242 }
243 }
244
245 if ( (SpinP_switch==3 && ( SO_switch==1 || Hub_U_switch==1 || 1<=Constraint_NCS_switch
246 || Zeeman_NCS_switch==1 || Zeeman_NCO_switch==1 )) && spin<=1 ){
247
248 for (m=0; m<ian; m++){
249 for (n=0; n<jan; n++){
250
251 iCDM[spin][Mc_AN][h_AN][m][n] = Mix_wgt1 * iCDM[spin][Mc_AN][h_AN][m][n]
252 + Mix_wgt2 * iPDM[spin][Mc_AN][h_AN][m][n];
253
254 iP2DM[spin][Mc_AN][h_AN][m][n] = iPDM[spin][Mc_AN][h_AN][m][n];
255 iPDM[spin][Mc_AN][h_AN][m][n] = iCDM[spin][Mc_AN][h_AN][m][n];
256 }
257 }
258 }
259
260 }
261 }
262
263 /*---------- added by TOYODA 18/JAN/2010 */
264 if (g_exx_switch) {
265 EXX_Simple_Mixing_DM(exx, Mix_wgt1,
266 exx_CDM[spin], exx_PDM[spin], exx_P2DM[spin]);
267 }
268 /*---------- until here */
269 }
270
271 }
272