1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <time.h>
5 #include "openmx_common.h"
6 #include "mpi.h"
7
8
9
Mixing_DM(int MD_iter,int SCF_iter,int SCF_iter0,int SucceedReadingDMfile,double *** ReRhok,double *** ImRhok,double ** Residual_ReRhok,double ** Residual_ImRhok,double * ReVk,double * ImVk,double * ReRhoAtomk,double * ImRhoAtomk)10 double Mixing_DM(int MD_iter,
11 int SCF_iter,
12 int SCF_iter0,
13 int SucceedReadingDMfile,
14 double ***ReRhok,
15 double ***ImRhok,
16 double **Residual_ReRhok,
17 double **Residual_ImRhok,
18 double *ReVk,
19 double *ImVk,
20 double *ReRhoAtomk,
21 double *ImRhoAtomk)
22 {
23 int pSCF_iter,NumMix,NumSlide;
24 int spin,ct_AN,wan1,TNO1,h_AN,Gh_AN,wan2;
25 int TNO2,i,j,ian,jan,m,n;
26 int spinmax,k1,k2,k3;
27 int Mc_AN,Gc_AN;
28 double time0;
29 double TStime,TEtime;
30 int numprocs,myid,ID;
31
32 /* MPI */
33 MPI_Comm_size(mpi_comm_level1,&numprocs);
34 MPI_Comm_rank(mpi_comm_level1,&myid);
35
36 MPI_Barrier(mpi_comm_level1);
37 dtime(&TStime);
38
39 /* In case of rmm-diish, then return. */
40 if (Cnt_switch!=1 && Mixing_switch==5) return 0.0;
41
42 /* In case of (rmm-diis || rmm-adiis) && SCF_iter<(Pulay_SCF/2)), then return. */
43 if ( Cnt_switch!=1 && SucceedReadingDMfile!=1
44 && (Mixing_switch==1 || Mixing_switch==6) && SCF_iter<=(Pulay_SCF/2) && Solver!=4) return 0.0;
45
46 /*******************************************************
47 Simple Mixing
48 *******************************************************/
49
50 if (Mixing_switch==0){
51
52 if (MD_iter==1 && SCF_iter==1){
53
54 /*---------- modified by TOYODA 18/JAN/2010 */
55
56 #if EXX_MIX_DM
57
58 if (g_exx_switch) {
59 Simple_Mixing_DM(0,1.0,DM[0],DM[1],DM[2],
60 iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2],
61 g_exx, g_exx_DM[0], g_exx_DM[1], g_exx_DM[2]);
62 } else {
63 Simple_Mixing_DM(0,1.0,DM[0],DM[1],DM[2],
64 iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2],
65 NULL, NULL, NULL, NULL);
66 }
67 #else
68 Simple_Mixing_DM( 0,1.0,DM[0],DM[1],DM[2],
69 iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2]);
70 #endif
71 /*---------- until here */
72
73 }
74 else{
75
76 /*---------- modified by TOYODA 18/JAN/2010 */
77 #if EXX_MIX_DM
78 if (g_exx_switch) {
79 Simple_Mixing_DM(1,Mixing_weight,DM[0],DM[1],DM[2],
80 iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2],
81 g_exx, g_exx_DM[0], g_exx_DM[1], g_exx_DM[2]);
82 } else {
83 Simple_Mixing_DM(1,Mixing_weight,DM[0],DM[1],DM[2],
84 iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2],
85 NULL, NULL, NULL, NULL);
86 }
87 #else
88 Simple_Mixing_DM(1,Mixing_weight,DM[0],DM[1],DM[2],iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2]);
89 #endif
90 /*---------- until here */
91
92 }
93 }
94
95 /*******************************************************
96 Pulay's method
97
98 Residual Minimazation Method (RMM) using
99 Direct Inversion in the Iterative Subspace (DIIS)
100 *******************************************************/
101
102 else if (Mixing_switch==1){
103
104 if (SCF_iter==1){
105 DIIS_Mixing_DM(1,ResidualDM,iResidualDM);
106 }
107 else if (SCF_iter<=(Pulay_SCF-1)){
108
109 /*---------- added by TOYODA 18/JAN/2010 */
110 #if EXX_MIX_DM
111 if (g_exx_switch) {
112 Simple_Mixing_DM(1,Mixing_weight,DM[0],DM[1],DM[2],
113 iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2],
114 g_exx, g_exx_DM[0], g_exx_DM[1], g_exx_DM[2]);
115 } else {
116 Simple_Mixing_DM(1,Mixing_weight,DM[0],DM[1],DM[2],
117 iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2],
118 NULL, NULL, NULL, NULL);
119 }
120 #else
121 Simple_Mixing_DM(1,Mixing_weight,DM[0],DM[1],DM[2],iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2]);
122 #endif
123 /*---------- until here */
124
125 }
126 else if (SCF_iter==Pulay_SCF)
127 DIIS_Mixing_DM(2,ResidualDM,iResidualDM);
128 else
129 DIIS_Mixing_DM(SCF_iter-(Pulay_SCF-2),ResidualDM,iResidualDM);
130 }
131
132 /*********************************************************
133 Guaranteed-Reduction Pulay's method (GR-Pulay method)
134 *********************************************************/
135
136 else if (Mixing_switch==2){
137 if (SCF_iter==1){
138 GR_Pulay_DM(1,ResidualDM);
139 }
140 else if (SCF_iter<=(Pulay_SCF-1))
141 /*---------- added by TOYODA 18/JAN/2010 */
142 #if EXX_MIX_DM
143 if (g_exx_switch) {
144 Simple_Mixing_DM(1,Mixing_weight,DM[0],DM[1],DM[2],
145 iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2],
146 g_exx, g_exx_DM[0], g_exx_DM[1], g_exx_DM[2]);
147 } else {
148 Simple_Mixing_DM(1,Mixing_weight,DM[0],DM[1],DM[2],
149 iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2],
150 NULL, NULL, NULL, NULL);
151 }
152 #else
153 Simple_Mixing_DM(1,Mixing_weight,DM[0],DM[1],DM[2],iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2]);
154 #endif
155 /*---------- until here */
156 else if (SCF_iter==Pulay_SCF)
157 GR_Pulay_DM(2,ResidualDM);
158 else
159 GR_Pulay_DM(SCF_iter-(Pulay_SCF-2),ResidualDM);
160 }
161
162 /*********************************************************
163 Kerker simple mixing in k-space
164 *********************************************************/
165
166 else if (Mixing_switch==3){
167
168 if (MD_iter==1 && SCF_iter0==1 && SucceedReadingDMfile==0){
169 Kerker_Mixing_Rhok(0, Mixing_weight,
170 ReRhok,ImRhok,
171 Residual_ReRhok,Residual_ImRhok,
172 ReVk,ImVk,ReRhoAtomk,ImRhoAtomk);
173 }
174
175 else if (MD_iter==1 && SCF_iter==1){
176 Kerker_Mixing_Rhok(0, Mixing_weight,
177 ReRhok,ImRhok,
178 Residual_ReRhok,Residual_ImRhok,
179 ReVk,ImVk,ReRhoAtomk,ImRhoAtomk);
180 }
181
182 else if ( (SCF_iter%30)==1) {
183 Kerker_Mixing_Rhok(1,Max_Mixing_weight,
184 ReRhok,ImRhok,
185 Residual_ReRhok,Residual_ImRhok,
186 ReVk,ImVk,ReRhoAtomk,ImRhoAtomk);
187 }
188
189 else{
190 Kerker_Mixing_Rhok(1,Mixing_weight,
191 ReRhok,ImRhok,
192 Residual_ReRhok,Residual_ImRhok,
193 ReVk,ImVk,ReRhoAtomk,ImRhoAtomk);
194 }
195 }
196
197 /*********************************************************
198 RMM-DIIS mixing with Kerker's weighting factor in k-space
199 *********************************************************/
200
201 else if (Mixing_switch==4){
202
203 if (MD_iter==1 && SCF_iter==1){
204 Kerker_Mixing_Rhok(0,Mixing_weight,
205 ReRhok,ImRhok,
206 Residual_ReRhok,Residual_ImRhok,
207 ReVk,ImVk,ReRhoAtomk,ImRhoAtomk);
208 }
209
210 else if (SCF_iter<Pulay_SCF){
211 Kerker_Mixing_Rhok(1,Mixing_weight,
212 ReRhok,ImRhok,
213 Residual_ReRhok,Residual_ImRhok,
214 ReVk,ImVk,ReRhoAtomk,ImRhoAtomk);
215 }
216
217 else{
218
219 if (6<SCF_iter){
220 DIIS_Mixing_Rhok(SCF_iter-(Pulay_SCF-3),
221 Mixing_weight,
222 ReRhok,ImRhok,
223 Residual_ReRhok,Residual_ImRhok,
224 ReVk,ImVk,ReRhoAtomk,ImRhoAtomk);
225 }
226 else{
227 DIIS_Mixing_Rhok(SCF_iter,
228 Mixing_weight,
229 ReRhok,ImRhok,
230 Residual_ReRhok,Residual_ImRhok,
231 ReVk,ImVk,ReRhoAtomk,ImRhoAtomk);
232 }
233
234 }
235 }
236
237 /*******************************************************
238 ADIIS method:
239 Augmented Roothaan-Hall (ARH) energy function
240 +
241 Direct Inversion in the Iterative Subspace (DIIS)
242 *******************************************************/
243
244 else if (Mixing_switch==6){
245
246 if (SCF_iter<=(Pulay_SCF-1)){
247 Simple_Mixing_DM(1,Mixing_weight,DM[0],DM[1],DM[2],
248 iDM[0],iDM[1],iDM[2],ResidualDM[2],iResidualDM[2],
249 NULL, NULL, NULL, NULL);
250 }
251 else
252 ADIIS_Mixing_DM(SCF_iter-Pulay_SCF+2,ResidualDM,iResidualDM);
253 }
254
255 /*******************************************************
256 not supported
257 *******************************************************/
258
259 else {
260 printf("Mixing_switch=%i is not supported.\n",Mixing_switch);
261 MPI_Finalize();
262 exit(0);
263 }
264
265 /* if SCF_iter0==1, then NormRD[0]=1 */
266
267 if (SCF_iter0==1) NormRD[0] = 1.0;
268
269 MPI_Barrier(mpi_comm_level1);
270 dtime(&TEtime);
271 time0 = TEtime - TStime;
272 return time0;
273 }
274
275
276
277