From b61b8e51e3eb57845a708a80adb138a6bf84b531 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 15 May 2009 01:01:27 +0200 Subject: [PATCH] Parallel and OK. --- Makefile | 4 +- debug.irp.f | 5 ++ eplf.irp.f | 24 +++---- eplf_grid.irp.f | 184 ++++++++++++++++++++++++++++++++++++++++++++++++ eplf_hf.irp.f | 5 ++ mpi.irp.f | 35 +++++++++ 6 files changed, 241 insertions(+), 16 deletions(-) create mode 100644 eplf_grid.irp.f create mode 100644 eplf_hf.irp.f create mode 100644 mpi.irp.f diff --git a/Makefile b/Makefile index 18b10b7..559ec11 100644 --- a/Makefile +++ b/Makefile @@ -1,5 +1,5 @@ -IRPF90 = irpf90 #-a -d -FC = ifort +IRPF90 = irpf90 -DMPI #-a -d +FC = mpif90 FCFLAGS= -O3 -xT SRC= diff --git a/debug.irp.f b/debug.irp.f index f2fa12a..ee3edb2 100644 --- a/debug.irp.f +++ b/debug.irp.f @@ -39,4 +39,9 @@ program debug print *, 'EPLF gamma : ', eplf_gamma print *, 'EPLF integral :', eplf_integral(i,j,eplf_gamma,point) print *, 'EPLF integral N :', eplf_integral_numeric(i,j,eplf_gamma,point) + + print *, '' + print *, 'EPLF grid Npoints :', grid_eplf_x_num, grid_eplf_y_num, grid_eplf_z_num + print *, 'EPLF grid step :', grid_eplf_step(:) + print *, 'EPLF grid origin :', grid_eplf_origin(:) end diff --git a/eplf.irp.f b/eplf.irp.f index 596738c..ae0d8fe 100644 --- a/eplf.irp.f +++ b/eplf.irp.f @@ -3,7 +3,7 @@ BEGIN_PROVIDER [ real, eplf_gamma ] BEGIN_DOC ! Value of the gaussian for the EPLF END_DOC - eplf_gamma = 10. + eplf_gamma = 1000. END_PROVIDER BEGIN_PROVIDER [ double precision, eplf_integral_matrix, (ao_num,ao_num) ] @@ -29,7 +29,9 @@ END_PROVIDER END_DOC integer :: i, j, k, l, m, n + double precision :: thr + thr = 1.d-12 / eplf_gamma eplf_up_up = 0. eplf_up_dn = 0. @@ -40,35 +42,29 @@ END_PROVIDER double precision :: ao_mn ao_mn = eplf_integral_matrix(m,n) - if (ao_mn /= 0.d0) then + if (abs(ao_mn) > thr) then do k=1,ao_num do l=1,ao_num double precision :: ao_klmn ao_klmn = ao_mn*ao_value_p(k)*ao_value_p(l) - if (ao_klmn /= 0.d0) then + if (abs(ao_klmn) > thr) then - do j=1,elec_beta_num + do j=1,elec_alpha_num if (mo_coef_transp(j,n) /= 0.d0) then - do i=1,elec_beta_num - eplf_up_up = eplf_up_up + 2.d0*ao_klmn* & + do i=1,elec_alpha_num + eplf_up_up = eplf_up_up + ao_klmn* & mo_coef_transp(i,k)*mo_coef_transp(j,n) * & (mo_coef_transp(i,l)*mo_coef_transp(j,m) - & mo_coef_transp(i,m)*mo_coef_transp(j,l) ) enddo endif enddo - - do j=1+elec_beta_num, elec_alpha_num + + do j=1,elec_beta_num if (mo_coef_transp(j,n) /= 0.d0) then do i=1,elec_beta_num - eplf_up_up = eplf_up_up + 2.d0*ao_klmn* & - mo_coef_transp(i,k)*mo_coef_transp(j,n) * & - (mo_coef_transp(i,l)*mo_coef_transp(j,m) - & - mo_coef_transp(i,m)*mo_coef_transp(j,l) ) - enddo - do i=1+elec_beta_num,elec_alpha_num eplf_up_up = eplf_up_up + ao_klmn* & mo_coef_transp(i,k)*mo_coef_transp(j,n) * & (mo_coef_transp(i,l)*mo_coef_transp(j,m) - & diff --git a/eplf_grid.irp.f b/eplf_grid.irp.f new file mode 100644 index 0000000..27173fe --- /dev/null +++ b/eplf_grid.irp.f @@ -0,0 +1,184 @@ +BEGIN_PROVIDER [ character*(100), grid_data_filename ] + BEGIN_DOC +! Name of the file containing the parameters of the grid + END_DOC + grid_data_filename = 'eplf_grid.data' +END_PROVIDER + +BEGIN_PROVIDER [ character*(100), grid_cube_filename ] + BEGIN_DOC +! Name of the file containing the parameters of the grid + END_DOC + grid_cube_filename = 'eplf_grid.cube' +END_PROVIDER + + BEGIN_PROVIDER [ integer, grid_eplf_x_num ] +&BEGIN_PROVIDER [ integer, grid_eplf_y_num ] +&BEGIN_PROVIDER [ integer, grid_eplf_z_num ] +&BEGIN_PROVIDER [ real , grid_eplf_step , (3) ] +&BEGIN_PROVIDER [ real , grid_eplf_origin, (3) ] + + real, parameter :: UNDEFINED=123456789.123456789 + + BEGIN_DOC +! Number of grid points in x, y, z directions + END_DOC + + integer :: Npoints(3) + real :: step_size(3) + real :: origin(3) + real :: opposite(3) + + namelist /eplf_grid/ Npoints, step_size, origin, opposite + + Npoints (:) = 80 + origin (:) = UNDEFINED + opposite (:) = UNDEFINED + step_size(:) = UNDEFINED + + logical :: do_next + inquire(file=grid_data_filename,exist=do_next) + if (do_next) then + open(99,file=grid_data_filename,action='READ',status='OLD') + read(99,eplf_grid,end=90,err=80) + close(99) + endif + + 90 continue + + if (origin(1) == UNDEFINED) then + integer :: i,l + do l=1,3 + origin(l) = nucl_coord(1,l) + do i=2,nucl_num + origin(l) = min(origin(l),nucl_coord(i,l)) + enddo + origin(l) = origin(l) - 4. + enddo + endif + + if (opposite(1) == UNDEFINED) then + do l=1,3 + opposite(l) = nucl_coord(1,l) + do i=2,nucl_num + opposite(l) = max(opposite(l),nucl_coord(i,l)) + enddo + opposite(l) = opposite(l) + 4. + enddo + endif + + if (step_size(1) == UNDEFINED) then + do l=1,3 + step_size(l) = (opposite(l) - origin(l))/float(Npoints(l)) + enddo + endif + + do l=1,3 + grid_eplf_origin(l) = origin(l) + grid_eplf_step(l) = step_size(l) + enddo + grid_eplf_x_num = Npoints(1) + grid_eplf_y_num = Npoints(2) + grid_eplf_z_num = Npoints(3) + + return + + 80 continue + print *, 'Error in input file' + stop 1 + +END_PROVIDER + +BEGIN_PROVIDER [ real, grid_eplf, (grid_eplf_x_num,grid_eplf_y_num,grid_eplf_z_num) ] + implicit none + BEGIN_DOC +! EPLF on a grid + END_DOC + + IRP_IF MPI + include 'mpif.h' + IRP_ENDIF + + integer :: ix, iy, iz + integer :: ibegin, iend + + do iz=1,grid_eplf_z_num + do iy=1,grid_eplf_y_num + do ix=1,grid_eplf_x_num + grid_eplf(ix,iy,iz) = 0. + enddo + enddo + enddo + + integer :: icount + icount = mpi_size + do iz=1,grid_eplf_z_num + if (mpi_master) then + print *, iz + endif + point(3) = grid_eplf_origin(3)+(iz-1)*grid_eplf_step(3) + do iy=1,grid_eplf_y_num + point(2) = grid_eplf_origin(2)+(iy-1)*grid_eplf_step(2) + do ix=1,grid_eplf_x_num + icount = icount-1 + if (icount == mpi_rank) then + point(1) = grid_eplf_origin(1)+(ix-1)*grid_eplf_step(1) + TOUCH point + grid_eplf(ix,iy,iz) = eplf_value + endif + if (icount == 0) then + icount = mpi_size + endif + enddo + enddo + enddo + + IRP_IF MPI + integer :: dim, ierr + real :: buffer(grid_eplf_x_num*grid_eplf_y_num*grid_eplf_z_num) + icount = 0 + do iz=1,grid_eplf_z_num + do iy=1,grid_eplf_y_num + do ix=1,grid_eplf_x_num + buffer(icount+ix) = grid_eplf(ix,iy,iz) + enddo + icount = icount + grid_eplf_x_num + enddo + enddo + dim = grid_eplf_x_num * grid_eplf_y_num * grid_eplf_z_num + call MPI_REDUCE(buffer,grid_eplf,dim,mpi_real, & + mpi_sum,0,MPI_COMM_WORLD,ierr) + call MPI_BARRIER(MPI_COMM_WORLD,ierr) + IRP_ENDIF + +END_PROVIDER + +subroutine write_grid_eplf + implicit none + integer :: i + integer :: l + integer :: ix, iy, iz + if (.not.mpi_master) then + return + endif + open(unit=99,file=grid_cube_filename,status='UNKNOWN',action='WRITE') + write (99,*) 'Cube File' + write (99,*) 'Analytical EPLF grid' + write (99,10) nucl_num,(grid_eplf_origin(i), i=1,3) + write (99,10) grid_eplf_x_num, grid_eplf_step(1), 0., 0. + write (99,10) grid_eplf_y_num, 0., grid_eplf_step(2), 0. + write (99,10) grid_eplf_z_num, 0., 0., grid_eplf_step(3) + do i=1,nucl_num + write (99,11) int(nucl_charge(i)), nucl_charge(i), (nucl_coord(i,l),l=1,3) + enddo + do ix = 1, grid_eplf_x_num + do iy = 1, grid_eplf_y_num + write (99,20) (grid_eplf(ix,iy,iz), iz=1, grid_eplf_z_num) + enddo + enddo + 10 format (2X,I3,3(2X,F10.6)) + 11 format (2X,I3,4(2X,F10.6)) + 20 format (6(E13.5)) + close(99) +end + diff --git a/eplf_hf.irp.f b/eplf_hf.irp.f new file mode 100644 index 0000000..b67765f --- /dev/null +++ b/eplf_hf.irp.f @@ -0,0 +1,5 @@ +program debug + PROVIDE ao_prim_num_max + call write_grid_eplf() +end + diff --git a/mpi.irp.f b/mpi.irp.f new file mode 100644 index 0000000..0f243cc --- /dev/null +++ b/mpi.irp.f @@ -0,0 +1,35 @@ + BEGIN_PROVIDER [ integer, mpi_rank ] +&BEGIN_PROVIDER [ integer, mpi_size ] +&BEGIN_PROVIDER [ logical, mpi_master ] + implicit none + + IRP_IF MPI + include 'mpif.h' + IRP_ENDIF + + BEGIN_DOC +! mpi_rank : ID of the current processor +! +! mpi_size : Total number of processors +! +! mpi_master : True if the current processor is the master + END_DOC + + + mpi_size = 1 + mpi_rank = 0 + + IRP_IF MPI + + integer :: ierr + call MPI_INIT(ierr) + call MPI_COMM_RANK(MPI_COMM_WORLD, mpi_rank, ierr) + call MPI_COMM_SIZE(MPI_COMM_WORLD, mpi_size, ierr) + + IRP_ENDIF + + mpi_master = (mpi_rank == 0) + +END_PROVIDER + +