PROGRAM wannier

IMPLICIT NONE 

CHARACTER(LEN=2),DIMENSION(:),ALLOCATABLE           :: element
REAL(KIND=8),DIMENSION(3)                           :: vec
INTEGER                                             :: stat                  
INTEGER                                             :: i, j, k, l, n, natom                                                         
REAL                                                :: dip_tot 
REAL(KIND=8),DIMENSION(:),ALLOCATABLE               :: dipole
REAL(KIND=8),DIMENSION(:,:),ALLOCATABLE             :: coord
    

OPEN(UNIT=50,FILE='h2o.xyz',STATUS='old',IOSTAT=stat) !reading number of atoms. You should save the optimized file as xyz format
READ(50,*) natom
CLOSE(50)

ALLOCATE(element(natom),coord(natom,3)) !alocating an array with the size of number of atoms
ALLOCATE(dipole(3))

dipole=0.

 OPEN(UNIT=51,FILE='h2o.xyz',STATUS='old',IOSTAT=stat) 
 READ(51,*) 
 READ(51,*)

 DO i=1,natom
 READ(51,*) element(i), coord(i,1), coord(i,2), coord(i,3) !reading elements and coordinates
 ENDDO

   DO i=1,natom
    IF(element(i).NE."O") CYCLE
     DO j=1,natom
      IF(element(j)=="O") CYCLE
       vec(:)=coord(j,:)-coord(i,:)
       IF(element(j) =="H") THEN
       dipole(:)=dipole(:) + vec*0.374710*1.60d-29 !charge*unit conversion
      ENDIF
     ENDDO
   ENDDO
    dip_tot=SQRT(dipole(1)**2+dipole(2)**2+dipole(3)**2)
  CLOSE(51)
  
  OPEN(UNIT=60,FILE='dipole_h2o.dat',STATUS='unknown',IOSTAT=stat)
  WRITE(60,*) REAL(dip_tot/3.33564d-30,KIND=8) !divided by unit conversion
  CLOSE(60)
    
DEALLOCATE(element,coord)
DEALLOCATE(dipole)

END PROGRAM wannier