fix some minor issues in TDDFT

git-svn-id: svn+ssh://chips.ncsu.edu/home/svn/repos/codes@3858 eb38b3b0-33ea-0310-94ad-d43fea0d214f

NEGF/current_scan.cpp

@@ -41,7 +41,7 @@ int main (int argc, char **argv)
     std::string whitespace_delims = " \n\t\"";
     std::vector<double> ener, cond;
     std::vector<std::string> line_values;
-    double bias = 0.05, temperature = 300.0;
+    double bias = 0.10, temperature = 300.0;
 
     double KT ;
     KT = temperature * Kb;

ON/Common/Main.cpp

@@ -131,6 +131,8 @@ using namespace El;
 #endif
 using namespace std;
 
+static void dipole_calculation(double *rhooo, double *dipole);
+
 int main(int argc, char **argv)
 {
 
@@ -255,6 +257,10 @@ int main(int argc, char **argv)
         RmgTimer *RTw = new RmgTimer("1-TOTAL: write");
 
         write_restart(ct.outfile, vh, vxc, vh_old, vxc_old, rho, rho_oppo, &states[0]); 
+        double dipole_ion, dipole_ele;
+        dipole_calculation(rhoc, &dipole_ion);
+        dipole_calculation(rho, &dipole_ele);
+        printf("\n duoile    %f %f \n", dipole_ion, dipole_ele);
 
         /* Save state information to file */
         // write_states_info(ct.outfile, &states[0]);
@@ -297,3 +303,45 @@ int main(int argc, char **argv)
     return 0;
 }
 
+
+static void dipole_calculation(double *rhooo, double *dipole)
+{
+
+
+
+    int i, j, k, idx;
+
+    double x, y, z;
+
+    dipole[0] = 0.0;
+    dipole[1] = 0.0;
+    dipole[2] = 0.0;
+
+    for(i = 0; i < get_FPX0_GRID(); i++)
+    {
+        x = (get_FPX_OFFSET() + i)*get_hxxgrid() * get_xside();
+        for(j = 0; j < get_FPY0_GRID(); j++)
+        {
+            y = (get_FPY_OFFSET() + j)*get_hyygrid() * get_yside();
+
+            for(k = 0; k < get_FPZ0_GRID(); k++)
+            {
+                z = (get_FPZ_OFFSET() + k)*get_hzzgrid() * get_zside();
+
+                idx = i * get_FPY0_GRID() * get_FPZ0_GRID() + j*get_FPZ0_GRID() + k;
+                dipole[0] += x * rhooo[idx];
+                dipole[1] += y * rhooo[idx];
+                dipole[2] += z * rhooo[idx];
+            }
+        }
+    }
+
+    dipole[0] *= get_vel_f();
+    dipole[1] *= get_vel_f();
+    dipole[2] *= get_vel_f();
+
+    idx = 3;
+    global_sums (dipole, &idx, pct.grid_comm);
+
+
+}

ON/Common/WriteRestart.cpp

@@ -103,10 +103,11 @@ void WriteRestart (char *name, double * vh, double *vxc, double *vh_old, double
 
 	/* write current ionic forces */
 	fprintf(fhandle,"\nionic_forces = \"");
+    for(int ic =0; ic < 4; ic++)
 	for (ion = 0; ion < ct.num_ions; ion++)
 	{
 	    iptr = &ct.ions[ion];
-	    fprintf(fhandle, "\n %#15.12g %#15.12g %#15.12g ", iptr->force[0][0], iptr->force[0][1], iptr->force[0][2]);
+	    fprintf(fhandle, "\n %#15.12g %#15.12g %#15.12g ", iptr->force[ic][0], iptr->force[ic][1], iptr->force[ic][2]);
 	}
 	fprintf(fhandle,"\n\"\n");
 

ON/Common/read_data.c

@@ -46,8 +46,6 @@ void read_data(char *name, double *vh, double *vxc, double *vh_old,
     my_barrier();
 
 
-    sprintf(newname, "%s%s", name, ".pot_rho");
-
     pe2xyz(pct.gridpe, &pex, &pey, &pez);
 
     int sizes[3], subsizes[3], starts[3];
@@ -84,19 +82,38 @@ void read_data(char *name, double *vh, double *vxc, double *vh_old,
 
     int amode = MPI_MODE_RDWR|MPI_MODE_CREATE;
     MPI_File mpi_fhand ;
+
+    sprintf(newname, "%s%s", name, ".vh");
     MPI_File_open(pct.grid_comm, newname, amode, fileinfo, &mpi_fhand);
+    disp=0;
+    MPI_File_set_view(mpi_fhand, disp, MPI_DOUBLE, filetype, "native", MPI_INFO_NULL);
+    MPI_File_read_all(mpi_fhand, vh, get_FP0_BASIS(),MPI_DOUBLE, &status);
+    MPI_File_close(&mpi_fhand);
+    MPI_Barrier(pct.grid_comm);
 
+    sprintf(newname, "%s%s", name, ".vxc");
+    MPI_File_open(pct.grid_comm, newname, amode, fileinfo, &mpi_fhand);
+    disp=0;
+    MPI_File_set_view(mpi_fhand, disp, MPI_DOUBLE, filetype, "native", MPI_INFO_NULL);
+    MPI_File_read_all(mpi_fhand, vxc, get_FP0_BASIS(),MPI_DOUBLE, &status);
+    MPI_File_close(&mpi_fhand);
+    MPI_Barrier(pct.grid_comm);
 
+    sprintf(newname, "%s%s", name, ".rho");
+    MPI_File_open(pct.grid_comm, newname, amode, fileinfo, &mpi_fhand);
     disp=0;
     MPI_File_set_view(mpi_fhand, disp, MPI_DOUBLE, filetype, "native", MPI_INFO_NULL);
+    MPI_File_read_all(mpi_fhand, rho, get_FP0_BASIS(),MPI_DOUBLE, &status);
+    MPI_File_close(&mpi_fhand);
+    MPI_Barrier(pct.grid_comm);
 
-    MPI_File_read(mpi_fhand, vh, get_FP0_BASIS(),MPI_DOUBLE, &status);
-    MPI_File_read(mpi_fhand, vxc, get_FP0_BASIS(),MPI_DOUBLE, &status);
-    MPI_File_read(mpi_fhand, rho, get_FP0_BASIS(),MPI_DOUBLE, &status);
-    MPI_File_read(mpi_fhand, vh_old, get_FP0_BASIS(),MPI_DOUBLE, &status);
-    MPI_File_read(mpi_fhand, vxc_old, get_FP0_BASIS(),MPI_DOUBLE, &status);
-    MPI_File_read(mpi_fhand, vh_corr, get_FP0_BASIS(),MPI_DOUBLE, &status);
+    sprintf(newname, "%s%s", name, ".vh_corr");
+    MPI_File_open(pct.grid_comm, newname, amode, fileinfo, &mpi_fhand);
+    disp=0;
+    MPI_File_set_view(mpi_fhand, disp, MPI_DOUBLE, filetype, "native", MPI_INFO_NULL);
+    MPI_File_read_all(mpi_fhand, vh_corr, get_FP0_BASIS(),MPI_DOUBLE, &status);
     MPI_File_close(&mpi_fhand);
+    MPI_Barrier(pct.grid_comm);
 
     my_barrier();
 

ON/Common/write_data.c

@@ -66,7 +66,6 @@ void write_data(char *name, double *vh, double *vxc, double *vh_old,
 	if(pct.gridpe == 0) rmg_printf("\n    Write start at %s\n", timeptr);
 	if(pct.gridpe == 0) fflush(NULL);
 
-	sprintf(newname, "%s%s", name, ".pot_rho");
 
 	pe2xyz (pct.gridpe, &pex, &pey, &pez);
 
@@ -104,19 +103,40 @@ void write_data(char *name, double *vh, double *vxc, double *vh_old,
 
 	amode = MPI_MODE_RDWR|MPI_MODE_CREATE;
 	MPI_File mpi_fhand ;
-	MPI_File_open(pct.grid_comm, newname, amode, fileinfo, &mpi_fhand);
-
 
+    MPI_Barrier(pct.grid_comm);
+	sprintf(newname, "%s%s", name, ".vh");
+	MPI_File_open(pct.grid_comm, newname, amode, fileinfo, &mpi_fhand);
 	disp=0;
 	MPI_File_set_view(mpi_fhand, disp, MPI_DOUBLE, filetype, "native", MPI_INFO_NULL);
-
 	MPI_File_write_all(mpi_fhand, vh, get_FP0_BASIS(),MPI_DOUBLE, &status);
+	MPI_File_close(&mpi_fhand);
+    MPI_Barrier(pct.grid_comm);
+
+	sprintf(newname, "%s%s", name, ".vxc");
+	MPI_File_open(pct.grid_comm, newname, amode, fileinfo, &mpi_fhand);
+	disp=0;
+	MPI_File_set_view(mpi_fhand, disp, MPI_DOUBLE, filetype, "native", MPI_INFO_NULL);
 	MPI_File_write_all(mpi_fhand, vxc, get_FP0_BASIS(),MPI_DOUBLE, &status);
+	MPI_File_close(&mpi_fhand);
+    MPI_Barrier(pct.grid_comm);
+
+	sprintf(newname, "%s%s", name, ".rho");
+	MPI_File_open(pct.grid_comm, newname, amode, fileinfo, &mpi_fhand);
+	disp=0;
+	MPI_File_set_view(mpi_fhand, disp, MPI_DOUBLE, filetype, "native", MPI_INFO_NULL);
 	MPI_File_write_all(mpi_fhand, rho, get_FP0_BASIS(),MPI_DOUBLE, &status);
-	MPI_File_write_all(mpi_fhand, vh_old, get_FP0_BASIS(),MPI_DOUBLE, &status);
-	MPI_File_write_all(mpi_fhand, vxc_old, get_FP0_BASIS(),MPI_DOUBLE, &status);
+	MPI_File_close(&mpi_fhand);
+    MPI_Barrier(pct.grid_comm);
+
+	sprintf(newname, "%s%s", name, ".vh_corr");
+	MPI_File_open(pct.grid_comm, newname, amode, fileinfo, &mpi_fhand);
+	disp=0;
+	MPI_File_set_view(mpi_fhand, disp, MPI_DOUBLE, filetype, "native", MPI_INFO_NULL);
 	MPI_File_write_all(mpi_fhand, vh_corr, get_FP0_BASIS(),MPI_DOUBLE, &status);
 	MPI_File_close(&mpi_fhand);
+    MPI_Barrier(pct.grid_comm);
+
 
 	my_barrier();
 

ON/Common/write_restart.c

@@ -101,6 +101,8 @@ void write_restart (char *name, double * vh, double *vxc, double *vh_old, double
 
 	/* write current ionic forces */
 	fprintf(fhandle,"\nionic_forces = \"");
+    int ic;
+    for (ic = 0; ic < 4; ic++)
 	for (ion = 0; ion < ct.num_ions; ion++)
 	{
 	    iptr = &ct.ions[ion];

ON/ON-NEGF-share/write_header.c

@@ -260,6 +260,9 @@ void write_header(void)
             case PSIPLOT:
                 printf("\n\n    PLOTTING PSI^2 IN DX FORM");
                 break;
+            case TDDFT:
+                printf("\n\n    TDDFT");
+                break;
 
             default:
                 error_handler("UNKNOWN MOLECULAR DYNAMICS METHOD");

TDDFT/Common/GetNewRho_rmgtddft.cpp

@@ -114,6 +114,7 @@ void GetNewRho_rmgtddft (double *psi, double *xpsi, double *rho, double *rho_mat
     rmg_printf ("normalization constant-1 for new charge is %f\n", t1-1);
     for(int i = 0;i < FP0_BASIS;i++) rho[i] *= t1;
 
+    delete [] rho_temp;
 
 
 }

TDDFT/Common/Main_rmg.cpp

@@ -357,9 +357,9 @@ template <typename OrbitalType> void run (Kpoint<OrbitalType> **Kptr)
         break;
 
     case TDDFT:
-       	Relax (0, vxc, vh, vnuc, rho, rho_oppo, rhocore, rhoc, Kptr);
+       //Relax (0, vxc, vh, vnuc, rho, rho_oppo, rhocore, rhoc, Kptr);
        // Quench (vxc, vh, vnuc, rho, rho_oppo, rhocore, rhoc, Kptr);
-        WriteRestart (ct.outfile, vh, rho, rho_oppo, vxc, Kptr);
+       // WriteRestart (ct.outfile, vh, rho, rho_oppo, vxc, Kptr);
         RmgTddft (vxc, vh, vnuc, rho, rho_oppo, rhocore, rhoc, Kptr);
         break;
     default:

TDDFT/Common/WriteData_rmgtddft.cpp

@@ -27,6 +27,7 @@
 #include <stdio.h>
 #if !(defined(_WIN32) || defined(_WIN64))
     #include <unistd.h>
+    #include <libgen.h>
 #else
     #include <io.h>
 #endif
@@ -46,7 +47,7 @@ void WriteData_rmgtddft (char *filename, double * vh, double * vxc,
 {
     int fhand;
     int fgrid_size;
-    char newname[MAX_PATH];
+    char newname[MAX_PATH], tmpname[MAX_PATH];
 
 
     int amode;
@@ -56,23 +57,46 @@ void WriteData_rmgtddft (char *filename, double * vh, double * vxc,
     fhand = open(newname, O_CREAT | O_TRUNC | O_RDWR, amode);
 
    if (fhand < 0) {
+
+       strcpy (tmpname, newname);
+
+
+#if !(defined(_WIN32) || defined(_WIN64))
+       if (!mkdir (dirname (tmpname), S_IRWXU))
+#else
+           char dirname[_MAX_DIR];
+       _splitpath(tmpname, NULL, dirname, NULL, NULL);
+       if (!_mkdir(dirname));
+#endif
+   if (1)
+        rmg_printf ("\n Creating directory %s succesfull\n\n", tmpname);
+    else
+        rmg_printf ("\n Creating directory %s FAILED\n\n", tmpname);
+
+
+    fhand = open(newname, O_CREAT | O_TRUNC | O_RDWR, amode);
+
+    if (fhand < 0) {
         rmg_printf("Can't open restart file %s", newname);
         rmg_error_handler(__FILE__, __LINE__, "Terminating.");
     }
 
 
-    fgrid_size = get_FPX0_GRID() * get_FPY0_GRID() * get_FPZ0_GRID();
-    write (fhand, vh, fgrid_size * sizeof(double));
-    write (fhand, vxc, fgrid_size * sizeof(double));
-    write (fhand, vh_corr, fgrid_size * sizeof(double));
+   }
+
+
+   fgrid_size = get_FPX0_GRID() * get_FPY0_GRID() * get_FPZ0_GRID();
+   write (fhand, vh, fgrid_size * sizeof(double));
+   write (fhand, vxc, fgrid_size * sizeof(double));
+   write (fhand, vh_corr, fgrid_size * sizeof(double));
+
+   int n2 = ct.num_states * ct.num_states;
 
-    int n2 = ct.num_states * ct.num_states;
-
-    write (fhand, Pn0, 2* n2 * sizeof(double));
-    write (fhand, Hmatrix, n2 * sizeof(double));
-    write (fhand, Smatrix, n2 * sizeof(double));
-    write (fhand, &tot_steps, sizeof(int));
-    close(fhand);
+   write (fhand, Pn0, 2* n2 * sizeof(double));
+   write (fhand, Hmatrix, n2 * sizeof(double));
+   write (fhand, Smatrix, n2 * sizeof(double));
+   write (fhand, &tot_steps, sizeof(int));
+   close(fhand);
 
 
 }                               /* end write_data */