scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

environ-related modifications, a few more input keywords for a new (under testing) environ feature. 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@10695 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
oliviero 2014-01-23 13:52:09 +00:00
parent 03a77bc34a
commit ad224920b4
4 changed files with 119 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
Modules/environ_input.f90
View File

@ -72,9 +72,9 @@ MODULE environ_input
! on atomic positions of width equal to atomicspread(ityp)
REAL(DP) :: solvationrad(nsx) = 3.D0
! solvationrad radius of the solvation shell for each species when the
! ionic dielectric function is adopted
! ionic dielectric function is adopted, in internal units (a.u.)
REAL(DP) :: atomicspread(nsx) = 0.5D0
! gaussian spreads of the atomic density of charge
! gaussian spreads of the atomic density of charge, in internal units (a.u.)
LOGICAL :: add_jellium = .false.
! depending on periodic boundary corrections, one may need to explicitly
! polarize the compensatinig jellium background
@ -123,6 +123,26 @@ MODULE environ_input
! density threshold for the onset of ionic countercharge
REAL(DP) :: solvent_temperature = 300.D0
! temperature of the solution
!
! External charges parameters
!
REAL(DP) :: env_extcharge_n = 0
! number of fixed external gaussian points/lines/planes of charges to be used
! in the calculation
REAL(DP) :: extcharge_origin(3) = 0.D0
! positions of the external charges are expressed with respect to this
! origin or (if no origin is specified) wrt the center of mass of the system
REAL(DP) :: extcharge_dim(nsx) = 0
! dimensionality of the external charge, 0=point, 1=line, 2=plane
REAL(DP) :: extcharge_axis(nsx) = 1
! axis along which the lines are placed or ortogonal to which the planes are placed
! 1=X, 2=Y, 3=Z
REAL(DP) :: extcharge_spread(nsx) = 0.5D0
! gaussian spread to be used in the generation of the ext. charges (in a.u.)
REAL(DP) :: extcharge_charge(nsx) = 0.D0
! charge of the external object in internal units (opposite of reality)
REAL(DP) :: extcharge_pos(3,nsx) = 0.D0
! position of the external object in internal units (a.u.)
NAMELIST / environ / &
verbose, environ_thr, environ_type, &
@ -134,7 +154,10 @@ MODULE environ_input
env_surface_tension, delta, &
env_pressure, &
env_ioncc_concentration, zion, rhopb, &
solvent_temperature
solvent_temperature, &
env_extcharge_n, extcharge_origin, extcharge_dim, &
extcharge_axis, extcharge_spread, extcharge_charge, &
extcharge_pos
CONTAINS
!
@ -186,6 +209,14 @@ MODULE environ_input
rhopb = 0.0001D0
solvent_temperature = 300.0D0
!
env_extcharge_n = 0
extcharge_origin(:) = 0.0D0
extcharge_dim(:) = 0
extcharge_axis(:) = 3
extcharge_spread(:) = 0.D0
extcharge_charge(:) = 0.D0
extcharge_pos(:,:) = 0.D0
!
RETURN
!
END SUBROUTINE
@ -239,6 +270,14 @@ MODULE environ_input
CALL mp_bcast( rhopb, ionode_id, intra_image_comm )
CALL mp_bcast( solvent_temperature, ionode_id, intra_image_comm )
!
CALL mp_bcast( env_extcharge_n, ionode_id, intra_image_comm )
CALL mp_bcast( extcharge_origin, ionode_id, intra_image_comm )
CALL mp_bcast( extcharge_dim, ionode_id, intra_image_comm )
CALL mp_bcast( extcharge_axis, ionode_id, intra_image_comm )
CALL mp_bcast( extcharge_charge, ionode_id, intra_image_comm )
CALL mp_bcast( extcharge_spread, ionode_id, intra_image_comm )
CALL mp_bcast( extcharge_pos, ionode_id, intra_image_comm )
!
RETURN
!
END SUBROUTINE

71
Modules/generate_function.f90
View File

@ -20,20 +20,22 @@ MODULE generate_function
CONTAINS
!----------------------------------------------------------------------
SUBROUTINE generate_gaussian( nnr, charge, spread, pos, rho )
SUBROUTINE generate_gaussian( nnr, dim, axis, charge, spread, pos, rho )
!----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE constants, ONLY : sqrtpi
USE cell_base, ONLY : at, bg, alat
USE io_global, ONLY : stdout
USE cell_base, ONLY : at, bg, alat, omega
USE fft_base, ONLY : dfftp
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : me_bgrp, intra_bgrp_comm
!
IMPLICIT NONE
!
! ... Declares variables
!
INTEGER, INTENT(IN) :: nnr
INTEGER, INTENT(IN) :: nnr, dim, axis
REAL( DP ), INTENT(IN) :: charge, spread
REAL( DP ), INTENT(IN) :: pos( 3 )
REAL( DP ), INTENT(INOUT) :: rho( nnr )
@ -44,7 +46,7 @@ CONTAINS
INTEGER :: index0
!
REAL( DP ) :: inv_nr1, inv_nr2, inv_nr3
REAL( DP ) :: scale, spr2, dist
REAL( DP ) :: scale, spr2, dist, lenght
REAL( DP ) :: r( 3 ), s( 3 )
REAL( DP ), ALLOCATABLE :: rholocal ( : )
!
@ -66,7 +68,19 @@ CONTAINS
ir_end = nnr
#endif
!
scale = charge / ( sqrtpi * spread )**3
IF (axis.LT.1.OR.axis.GT.3) &
WRITE(stdout,*)'WARNING: wrong axis in generate_gaussian'
IF ( dim .EQ. 0 ) THEN
scale = charge / ( sqrtpi * spread )**3
ELSE IF ( dim .EQ. 1 ) THEN
lenght = at(axis,axis) * alat
scale = charge / lenght / ( sqrtpi * spread )**2
ELSE IF ( dim .EQ. 2 ) THEN
lenght = at(axis,axis) * alat
scale = charge * lenght / omega / ( sqrtpi * spread )
ELSE
WRITE(stdout,*)'WARNING: wrong dim in generate_gaussian'
ENDIF
spr2 = ( spread / alat )**2
ALLOCATE( rholocal( nnr ) )
rholocal = 0.D0
@ -89,6 +103,16 @@ CONTAINS
!
r(:) = pos(:) - r(:)
!
! ... possibly 2D or 1D gaussians
!
IF ( dim .EQ. 1) THEN
r(axis) = 0.D0
ELSE IF ( dim .EQ. 2 ) THEN
DO i = 1, 3
IF ( i .NE. axis ) r(i) = 0.D0
ENDDO
END IF
!
! ... minimum image convention
!
s(:) = MATMUL( r(:), bg(:,:) )
@ -110,12 +134,13 @@ CONTAINS
END SUBROUTINE generate_gaussian
!----------------------------------------------------------------------
!----------------------------------------------------------------------
SUBROUTINE generate_gradgaussian( nnr, charge, spread, pos, gradrho )
SUBROUTINE generate_gradgaussian( nnr, dim, axis, charge, spread, pos, gradrho )
!----------------------------------------------------------------------
!
USE kinds, ONLY: DP
USE constants, ONLY: sqrtpi
USE cell_base, ONLY : at, bg, alat
USE kinds, ONLY : DP
USE constants, ONLY : sqrtpi
USE io_global, ONLY : stdout
USE cell_base, ONLY : at, bg, alat, omega
USE fft_base, ONLY : dfftp
USE mp_bands, ONLY : me_bgrp, intra_bgrp_comm
!
@ -123,7 +148,7 @@ CONTAINS
!
! ... Declares variables
!
INTEGER, INTENT(IN) :: nnr
INTEGER, INTENT(IN) :: nnr, dim, axis
REAL( DP ), INTENT(IN) :: charge, spread
REAL( DP ), INTENT(IN) :: pos( 3 )
REAL( DP ), INTENT(INOUT) :: gradrho( 3, nnr )
@ -134,7 +159,7 @@ CONTAINS
INTEGER :: index0
!
REAL( DP ) :: inv_nr1, inv_nr2, inv_nr3
REAL( DP ) :: scale, spr2, dist
REAL( DP ) :: scale, spr2, dist, lenght
REAL( DP ) :: r( 3 ), s( 3 )
REAL( DP ), ALLOCATABLE :: gradrholocal ( :, : )
!
@ -156,7 +181,19 @@ CONTAINS
ir_end = nnr
#endif
!
scale = 2.d0 * charge / sqrtpi**3 / spread**5
IF (axis.LT.1.OR.axis.GT.3) &
WRITE(stdout,*)'WARNING: wrong axis in generate_gaussian'
IF ( dim .EQ. 0 ) THEN
scale = charge / ( sqrtpi * spread )**3
ELSE IF ( dim .EQ. 1 ) THEN
lenght = at(axis,axis) * alat
scale = charge / lenght / ( sqrtpi * spread )**2
ELSE IF ( dim .EQ. 2 ) THEN
lenght = at(axis,axis) * alat
scale = charge * lenght / omega / ( sqrtpi * spread )
ELSE
WRITE(stdout,*)'WARNING: wrong dim in generate_gaussian'
ENDIF
spr2 = ( spread / alat )**2
ALLOCATE( gradrholocal( 3, nnr ) )
gradrholocal = 0.D0
@ -179,6 +216,16 @@ CONTAINS
!
r(:) = pos(:) - r(:)
!
! ... possibly 2D or 1D gaussians
!
IF ( dim .EQ. 1) THEN
r(axis) = 0.D0
ELSE IF ( dim .EQ. 2 ) THEN
DO i = 1, 3
IF ( i .NE. axis ) r(i) = 0.D0
ENDDO
END IF
!
! ... minimum image convention
!
s(:) = MATMUL( r(:), bg(:,:) )

11
PW/src/electrons.f90
View File

@ -318,8 +318,9 @@ SUBROUTINE electrons_scf ( no_printout )
env_static_permittivity, &
env_surface_tension, env_pressure, &
env_periodicity, env_ioncc_concentration,&
deenviron, esolvent, ecavity, epressure, &
eperiodic, eioncc
env_extcharge_n, deenviron, esolvent, &
ecavity, epressure, eperiodic, eioncc, &
eextcharge
#endif
USE dfunct, ONLY : newd
USE esm, ONLY : do_comp_esm, esm_printpot
@ -639,7 +640,7 @@ SUBROUTINE electrons_scf ( no_printout )
vltot = vltot_zero
!
CALL calc_eenviron( dfftp%nnr, nspin, rhoin%of_r, deenviron, esolvent, &
ecavity, epressure, eperiodic, eioncc )
ecavity, epressure, eperiodic, eioncc, eextcharge )
!
update_venviron = .NOT. conv_elec .AND. dr2 .LT. environ_thr
!
@ -733,7 +734,7 @@ SUBROUTINE electrons_scf ( no_printout )
! ... adds the external environment contribution to the energy
!
IF ( do_environ ) etot = etot + deenviron + esolvent + ecavity + &
epressure + eperiodic + eioncc
epressure + eperiodic + eioncc + eextcharge
#endif
!
IF ( .NOT. no_printout ) CALL print_energies ( )
@ -1109,6 +1110,7 @@ SUBROUTINE electrons_scf ( no_printout )
ELSE IF ( env_periodicity .NE. 3 ) THEN
WRITE( stdout, 9204 ) eperiodic
ENDIF
IF ( env_extcharge_n .GT. 0 ) WRITE( stdout, 9206 ) eextcharge
ENDIF
!
#endif
@ -1165,6 +1167,7 @@ SUBROUTINE electrons_scf ( no_printout )
9203 FORMAT( ' PV energy =',F17.8,' Ry' )
9204 FORMAT( ' periodic energy correct. =',F17.8,' Ry' )
9205 FORMAT( ' ionic charge energy =',F17.8,' Ry' )
9206 FORMAT( ' external charges energy =',F17.8,' Ry' )
#endif
END SUBROUTINE print_energies

14
PW/src/input.f90
View File

@ -241,7 +241,7 @@ SUBROUTINE iosys()
!
! ... ENVIRON namelist
!
USE environ_input, ONLY : verbose, environ_thr, environ_type, &
USE environ_input, ONLY : verbose, environ_thr, environ_type, &
stype, rhomax, rhomin, tbeta, &
env_static_permittivity, eps_mode, &
solvationrad, atomicspread, add_jellium, &
@ -250,7 +250,11 @@ SUBROUTINE iosys()
env_surface_tension, delta, &
env_pressure, &
env_ioncc_concentration, zion, rhopb, &
solvent_temperature
solvent_temperature, &
env_extcharge_n, extcharge_origin, &
extcharge_dim, extcharge_axis, &
extcharge_pos, extcharge_spread, &
extcharge_charge
#endif
!
! ... ELECTRONS namelist
@ -1285,7 +1289,11 @@ SUBROUTINE iosys()
env_surface_tension, delta, &
env_pressure, &
env_ioncc_concentration, zion, rhopb, &
solvent_temperature )
solvent_temperature, &
env_extcharge_n, extcharge_origin, &
extcharge_dim, extcharge_axis, &
extcharge_pos, extcharge_spread, &
extcharge_charge )
!
IF ( do_environ ) CALL environ_initions_allocate( nat_, ntyp )
!