MODULE neutralize
!
!**********************************************************************************
!*  NEUTRALIZE                                                                    *
!**********************************************************************************
!* This module verifies if the system is "neutral", i.e. that the sum of the      *
!* given property is zero. For instance, if the given property is electric charge *
!* then this module verifies if the system is charge-neutral.                     *
!* If it is not, then it removes atoms until it becomes neutral (or close).       *
!**********************************************************************************
!* (C) April 2015 - Pierre Hirel                                                  *
!*     Unité Matériaux Et Transformations (UMET),                                 *
!*     Université de Lille 1, Bâtiment C6, F-59655 Villeneuve D'Ascq (FRANCE)     *
!*     pierre.hirel@univ-lille.fr                                                 *
!* Last modification: P. Hirel - 08 April 2015                                    *
!**********************************************************************************
!* This program is free software: you can redistribute it and/or modify           *
!* it under the terms of the GNU General Public License as published by           *
!* the Free Software Foundation, either version 3 of the License, or              *
!* (at your option) any later version.                                            *
!*                                                                                *
!* This program is distributed in the hope that it will be useful,                *
!* but WITHOUT ANY WARRANTY; without even the implied warranty of                 *
!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                  *
!* GNU General Public License for more details.                                   *
!*                                                                                *
!* You should have received a copy of the GNU General Public License              *
!* along with this program.  If not, see <http://www.gnu.org/licenses/>.          *
!**********************************************************************************
!
USE comv
USE constants
USE messages
USE neighbors
USE files
USE subroutines
!
!
CONTAINS
!
!
SUBROUTINE NEUTRALIZE_XYZ(H,P,S,AUXNAMES,AUX,SELECT,property)
!
!
IMPLICIT NONE
CHARACTER(LEN=2):: species
CHARACTER(LEN=128):: msg
CHARACTER(LEN=128),INTENT(IN):: property
CHARACTER(LEN=128),DIMENSION(:),ALLOCATABLE,INTENT(INOUT):: AUXNAMES !names of auxiliary properties
LOGICAL:: exceeds100 !does the number of neighbours exceed 100?
LOGICAL,DIMENSION(:),ALLOCATABLE,INTENT(IN):: SELECT  !mask for atom list
INTEGER:: i, j
INTEGER:: iprop
INTEGER:: Nremoved
INTEGER,DIMENSION(:),ALLOCATABLE:: Nlist !index of neighbours
REAL(dp):: SumProp  !sum of values for given property
REAL(dp),DIMENSION(3,3),INTENT(IN):: H   !Base vectors of the supercell
REAL(dp),DIMENSION(:,:),ALLOCATABLE:: V_NN  !final positions of 1st NN
REAL(dp),DIMENSION(:,:),ALLOCATABLE,INTENT(INOUT):: P, S !positions of atoms, shells
REAL(dp),DIMENSION(:,:),ALLOCATABLE:: Q, T               !positions of atoms, shells (temporary)
REAL(dp),DIMENSION(:,:),ALLOCATABLE,INTENT(INOUT):: AUX  !auxiliary properties of atoms
REAL(dp),DIMENSION(:,:),ALLOCATABLE:: newAUX             !auxiliary properties of atoms (temporary)
!
!
!Initialize variables
i = 0
iprop = 0
SumProp = 0.d0
IF(ALLOCATED(Q)) DEALLOCATE(Q)
IF(ALLOCATED(T)) DEALLOCATE(T)
IF(ALLOCATED(newAUX)) DEALLOCATE(newAUX)
!
!Locate the index of the given property
DO i=1,SIZE(AUXNAMES)
  IF( ADJUSTL(AUXNAMES(i)) == ADJUSTL(property) ) THEN
    iprop = i
  ENDIF
ENDDO
!
IF( iprop<=0 ) THEN
  !this property doesn't exist
ENDIF
!
!Compute the sum of the given property
SumProp = SUM(AUX(:,iprop))
!
!
!
100 CONTINUE
Nremoved = 0
IF( SumProp .NE. 0.d0 ) THEN
  !Find a way to make the system neutral
ENDIF
!
!
!
300 CONTINUE
IF( Nremoved==1 ) THEN
  CALL ATOMSK_MSG(2080,(/species/),(/DBLE(Nremoved),DBLE(SIZE(P,1))/))
ELSE
  CALL ATOMSK_MSG(2080,(/''/),(/DBLE(Nremoved),DBLE(SIZE(P,1))/))
ENDIF
GOTO 1000
!
!
!
800 CONTINUE
CALL ATOMSK_MSG(802,(/''/),(/DBLE(i)/))
nerr = nerr+1
!
!
!
1000 CONTINUE
!
!
END SUBROUTINE NEUTRALIZE_XYZ
!
!
END MODULE neutralize
