1 /**********************************************************************
2 Set_Initial_DM.c:
3
4 Set_Initial_DM.c is a subroutine to set an initial density matrix.
5
6 Log of Set_Initial_DM.c:
7
8 7/Feb./2012 Released by T.Ozaki
9
10 ***********************************************************************/
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <math.h>
16 #include <time.h>
17 #include "openmx_common.h"
18 #include "mpi.h"
19 #include <omp.h>
20
21 #define measure_time 0
22
23
Set_Initial_DM(double ***** CDM,double ***** H)24 double Set_Initial_DM(double *****CDM, double *****H)
25 {
26 int i,j,spin,spinmax,po,loopN;
27 int Mc_AN,Gc_AN,h_AN,Gh_AN,Cwan,Hwan;
28 double cp,cp_min,cp_max,ns,tns,dn,ff,x,TZ;
29 double spin_degeneracy;
30 double TStime,TEtime,time0;
31 int myid,numprocs;
32
33 MPI_Comm_size(mpi_comm_level1,&numprocs);
34 MPI_Comm_rank(mpi_comm_level1,&myid);
35
36 /* measure TStime */
37 MPI_Barrier(mpi_comm_level1);
38 dtime(&TStime);
39
40 /* initialize CDM */
41
42 for (spin=0; spin<=SpinP_switch; spin++){
43 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
44 Gc_AN = M2G[Mc_AN];
45 Cwan = WhatSpecies[Gc_AN];
46 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
47 Gh_AN = natn[Gc_AN][h_AN];
48 Hwan = WhatSpecies[Gh_AN];
49 for (i=0; i<Spe_Total_NO[Cwan]; i++){
50 for (j=0; j<Spe_Total_NO[Hwan]; j++){
51 CDM[spin][Mc_AN][h_AN][i][j] = 0.0;
52 }
53 }
54 }
55 }
56 }
57
58 /* set diagonal terms of CDM */
59
60 if (SpinP_switch==0){
61 spinmax = 0;
62 spin_degeneracy = 2.0;
63 }
64 else if (SpinP_switch==1){
65 spinmax = 1;
66 spin_degeneracy = 1.0;
67 }
68 else if (SpinP_switch==3){
69 spinmax = 1;
70 spin_degeneracy = 1.0;
71 }
72
73 TZ = 0.0;
74 for (i=1; i<=atomnum; i++){
75 Cwan = WhatSpecies[i];
76 TZ = TZ + Spe_Core_Charge[Cwan];
77 }
78
79 /* CDM[0] and CDM[1] */
80
81 cp_max = 30.0;
82 cp_min = -30.0;
83 loopN = 0;
84 po = 0;
85
86 do {
87
88 cp = 0.50*(cp_max + cp_min);
89 ns = 0.0;
90
91 for (spin=0; spin<=spinmax; spin++){
92 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
93
94 Gc_AN = M2G[Mc_AN];
95 Cwan = WhatSpecies[Gc_AN];
96
97 for (i=0; i<Spe_Total_NO[Cwan]; i++){
98 x = (H[spin][Mc_AN][0][i][i] - cp)*Beta;
99 ff = 1.0/(1.0 + exp(x));
100 CDM[spin][Mc_AN][0][i][i] = ff;
101 ns += spin_degeneracy*ff;
102 }
103
104 /*
105 for (i=0; i<Spe_Total_NO[Cwan]; i++){
106 printf("Gc_AN=%2d i=%2d CDM=%15.12f\n",Gc_AN,i,CDM[spin][Mc_AN][0][i][i]);
107 }
108 */
109
110 } /* Mc_AN */
111 } /* spin */
112
113 MPI_Allreduce(&ns, &tns, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
114
115 dn = TZ - tns - system_charge;
116
117 if (0.0<=dn) cp_min = cp;
118 else cp_max = cp;
119
120 if (fabs(dn)<1.0e-12) po = 1;
121
122 printf("loopN=%3d cp=%15.12f TZ=%15.12f dn=%15.12f\n",loopN,cp,TZ,dn);
123
124 loopN++;
125
126 } while (po==0 && loopN<1000);
127
128 /* CDM[2] */
129 if (SpinP_switch==3){
130 }
131
132 /* calculate time0 */
133 MPI_Barrier(mpi_comm_level1);
134 dtime(&TEtime);
135 time0 = TEtime - TStime;
136
137 return time0;
138 }
139