Commit 4cd5b32e authored by oliviero's avatar oliviero

Modified a bit the Environ structure and added a further contribution...

Modified a bit the Environ structure and added a further contribution (correction of slab boundary conditions).
All modifications within #ifdef __ENVIRON #endif statements, no effects on actual code.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@8702 c92efa57-630b-4861-b058-cf58834340f0
parent 0c7197e2
......@@ -5,32 +5,29 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!
!----------------------------------------------------------------------
! Module to generate functions on the real space dense grid
! Written by Oliviero Andreussi
!----------------------------------------------------------------------
!
!=----------------------------------------------------------------------=!
MODULE generate_function
MODULE generate_function
!=----------------------------------------------------------------------=!
USE kinds, ONLY: DP
USE kinds, ONLY: DP
IMPLICIT NONE
IMPLICIT NONE
CONTAINS
!----------------------------------------------------------------------------
CONTAINS
!----------------------------------------------------------------------
SUBROUTINE generate_gaussian( nrxx, charge, spread, pos, rho )
!----------------------------------------------------------------------------
!----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE constants, ONLY : sqrtpi
USE cell_base, ONLY : at, bg, alat
USE fft_base, ONLY : dfftp
USE mp_global, ONLY : me_pool, intra_pool_comm, &
me_bgrp, intra_bgrp_comm
USE mp_global, ONLY : me_bgrp, intra_bgrp_comm
!
IMPLICIT NONE
!
......@@ -58,17 +55,9 @@
index0 = 0
!
#if defined (__MPI)
!
#ifdef __BANDS
DO i = 1, me_bgrp
index0 = index0 + dfftp%nr1x*dfftp%nr2x*dfftp%npp(i)
END DO
#else
DO i = 1, me_pool
index0 = index0 + dfftp%nr1x*dfftp%nr2x*dfftp%npp(i)
END DO
#endif
!
#endif
!
scale = charge / ( sqrtpi * spread )**3
......@@ -112,19 +101,18 @@
!
RETURN
!
!----------------------------------------------------------------------------
!----------------------------------------------------------------------
END SUBROUTINE generate_gaussian
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
!----------------------------------------------------------------------
!----------------------------------------------------------------------
SUBROUTINE generate_gradgaussian( nrxx, charge, spread, pos, gradrho )
!----------------------------------------------------------------------------
!----------------------------------------------------------------------
!
USE kinds, ONLY: DP
USE constants, ONLY: sqrtpi
USE cell_base, ONLY : at, bg, alat
USE fft_base, ONLY : dfftp
USE mp_global, ONLY : me_pool, intra_pool_comm, &
me_bgrp, intra_bgrp_comm
USE mp_global, ONLY : me_bgrp, intra_bgrp_comm
!
IMPLICIT NONE
!
......@@ -152,17 +140,9 @@
index0 = 0
!
#if defined (__MPI)
!
#ifdef __BANDS
DO i = 1, me_bgrp
index0 = index0 + dfftp%nr1x*dfftp%nr2x*dfftp%npp(i)
END DO
#else
DO i = 1, me_pool
index0 = index0 + dfftp%nr1x*dfftp%nr2x*dfftp%npp(i)
END DO
#endif
!
#endif
!
scale = 2.d0 * charge / sqrtpi**3 / spread**5
......@@ -206,18 +186,17 @@
!
RETURN
!
!--------------------------------------------------------------------
!----------------------------------------------------------------------
END SUBROUTINE generate_gradgaussian
!--------------------------------------------------------------------
!----------------------------------------------------------------------------
!----------------------------------------------------------------------
!----------------------------------------------------------------------
SUBROUTINE generate_exponential( nrxx, spread, pos, rho )
!----------------------------------------------------------------------------
!----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg, alat
USE fft_base, ONLY : dfftp
USE mp_global, ONLY : me_pool, intra_pool_comm, &
me_bgrp, intra_bgrp_comm
USE mp_global, ONLY : me_bgrp, intra_bgrp_comm
!
IMPLICIT NONE
!
......@@ -246,17 +225,9 @@
index0 = 0
!
#if defined (__MPI)
!
#ifdef __BANDS
DO i = 1, me_bgrp
index0 = index0 + dfftp%nr1x*dfftp%nr2x*dfftp%npp(i)
END DO
#else
DO i = 1, me_pool
index0 = index0 + dfftp%nr1x*dfftp%nr2x*dfftp%npp(i)
END DO
#endif
!
#endif
!
ALLOCATE( rholocal( nrxx ) )
......@@ -304,18 +275,17 @@
!
RETURN
!
!--------------------------------------------------------------------
!----------------------------------------------------------------------
END SUBROUTINE generate_exponential
!--------------------------------------------------------------------
!----------------------------------------------------------------------------
!----------------------------------------------------------------------
!----------------------------------------------------------------------
SUBROUTINE generate_gradexponential( nrxx, spread, pos, gradrho )
!----------------------------------------------------------------------------
!----------------------------------------------------------------------
!
USE kinds, ONLY: DP
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg, alat
USE fft_base, ONLY : dfftp
USE mp_global, ONLY : me_pool, intra_pool_comm, &
me_bgrp, intra_bgrp_comm
USE mp_global, ONLY : me_bgrp, intra_bgrp_comm
!
IMPLICIT NONE
!
......@@ -344,17 +314,9 @@
index0 = 0
!
#if defined (__MPI)
!
#ifdef __BANDS
DO i = 1, me_bgrp
index0 = index0 + dfftp%nr1x*dfftp%nr2x*dfftp%npp(i)
END DO
#else
DO i = 1, me_pool
index0 = index0 + dfftp%nr1x*dfftp%nr2x*dfftp%npp(i)
END DO
#endif
!
#endif
!
ALLOCATE( gradrholocal( 3, nrxx ) )
......@@ -400,11 +362,76 @@
!
RETURN
!
!--------------------------------------------------------------------
!----------------------------------------------------------------------
END SUBROUTINE generate_gradexponential
!--------------------------------------------------------------------
!----------------------------------------------------------------------
!----------------------------------------------------------------------
SUBROUTINE generate_axis( nrxx, icor, pos, axis )
!----------------------------------------------------------------------
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg, alat
USE fft_base, ONLY : dfftp
USE mp_global, ONLY : me_bgrp, intra_bgrp_comm
!
INTEGER, INTENT(IN) :: nrxx
INTEGER, INTENT(IN) :: icor
REAL(DP), INTENT(IN) :: pos(3)
REAL(DP), INTENT(OUT) :: axis( dfftp%nnr )
!
INTEGER :: i, j, k, ir, ip, index, index0
REAL(DP) :: inv_nr1, inv_nr2, inv_nr3
REAL(DP) :: r(3), s(3)
!
inv_nr1 = 1.D0 / DBLE( dfftp%nr1 )
inv_nr2 = 1.D0 / DBLE( dfftp%nr2 )
inv_nr3 = 1.D0 / DBLE( dfftp%nr3 )
!
index0 = 0
!
#if defined (__MPI)
DO i = 1, me_bgrp
index0 = index0 + dfftp%nr1x*dfftp%nr2x*dfftp%npp(i)
END DO
#endif
!
DO ir = 1, dfftp%nnr
!
! ... three dimensional indexes
!
index = index0 + ir - 1
k = index / (dfftp%nr1x*dfftp%nr2x)
index = index - (dfftp%nr1x*dfftp%nr2x)*k
j = index / dfftp%nr1x
index = index - dfftp%nr1x*j
i = index
!
DO ip = 1, 3
r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &
DBLE( j )*inv_nr2*at(ip,2) + &
DBLE( k )*inv_nr3*at(ip,3)
END DO
!
r(:) = r(:) - pos(:)
!
! ... minimum image convention
!
s(:) = MATMUL( r(:), bg(:,:) )
s(:) = s(:) - ANINT(s(:))
r(:) = MATMUL( at(:,:), s(:) )
!
axis(ir) = r(icor)
!
END DO
!
axis = axis * alat
!
RETURN
!
!----------------------------------------------------------------------
END SUBROUTINE generate_axis
!----------------------------------------------------------------------
!=----------------------------------------------------------------------=!
END MODULE generate_function
END MODULE generate_function
!=----------------------------------------------------------------------=!
......@@ -76,7 +76,8 @@ SUBROUTINE electrons()
vltot_zero, environ_thr, &
env_static_permittivity, &
env_surface_tension, env_pressure, &
deenviron, esolvent, ecavity, epressure
env_slab_geometry, deenviron, &
esolvent, ecavity, epressure, eslab
#endif
USE dfunct, only : newd
USE esm, ONLY : do_comp_esm, esm_printpot
......@@ -180,8 +181,9 @@ SUBROUTINE electrons()
#ifdef __ENVIRON
IF ( do_environ ) THEN
vltot_zero = vltot
CALL environ_initions( nat, nsp, ityp, zv, tau )
CALL environ_initcell( dfftp%nr1*dfftp%nr2*dfftp%nr3, omega )
CALL environ_initions( dfftp%nnr, nat, nsp, ityp, zv, tau )
CALL environ_initcell( dfftp%nnr, dfftp%nr1*dfftp%nr2*dfftp%nr3, &
omega, alat )
END IF
#endif
!
......@@ -455,7 +457,7 @@ SUBROUTINE electrons()
vltot = vltot_zero
!
CALL calc_eenviron( dfftp%nnr, nspin, rhoin%of_r, vltot_zero, &
deenviron, esolvent, ecavity, epressure )
deenviron, esolvent, ecavity, epressure, eslab )
!
update_venviron = .NOT. conv_elec .AND. dr2 .LT. environ_thr
!
......@@ -592,7 +594,7 @@ SUBROUTINE electrons()
!
! ... adds the external environment contribution to the energy
!
IF ( do_environ ) etot = etot + deenviron + esolvent + ecavity + epressure
IF ( do_environ ) etot = etot + deenviron + esolvent + ecavity + epressure + eslab
#endif
!
IF ( ( conv_elec .OR. MOD( iter, iprint ) == 0 ) .AND. .NOT. lmd ) THEN
......@@ -648,7 +650,9 @@ SUBROUTINE electrons()
IF ( env_static_permittivity .GT. 1.D0 ) WRITE( stdout, 9201 ) esolvent
IF ( env_surface_tension .GT. 0.D0 ) WRITE( stdout, 9202 ) ecavity
IF ( env_pressure .NE. 0.D0 ) WRITE( stdout, 9203 ) epressure
IF ( env_slab_geometry ) WRITE( stdout, 9204 ) eslab
ENDIF
!
#endif
!
IF ( lsda ) WRITE( stdout, 9017 ) magtot, absmag
......@@ -773,6 +777,7 @@ SUBROUTINE electrons()
9201 FORMAT( ' solvation energy =',F17.8,' Ry' )
9202 FORMAT( ' cavitation energy =',F17.8,' Ry' )
9203 FORMAT( ' PV energy =',F17.8,' Ry' )
9204 FORMAT( ' slab energy correction =',F17.8,' Ry' )
#endif
!
CONTAINS
......
......@@ -40,7 +40,7 @@ SUBROUTINE forces()
USE control_flags, ONLY : gamma_only, remove_rigid_rot, textfor, &
iverbosity, llondon
#ifdef __ENVIRON
USE environ_base, ONLY : do_environ, env_static_permittivity, eps_mode, rhopol
USE environ_base, ONLY : do_environ, env_static_permittivity, rhopol
#endif
USE bp, ONLY : lelfield, gdir, l3dstring, efield_cart, &
efield_cry,efield
......@@ -131,29 +131,23 @@ SUBROUTINE forces()
!
! ... The external environment contribution
!
CALL start_clock('calc_fsolv')
!
ALLOCATE ( force_environ ( 3 , nat ) )
!
force_environ = 0.0_DP
!
! ... Computes here the solvent contribution
!
IF ( env_static_permittivity .GT. 1.D0 ) &
CALL force_lc( nat, tau, ityp, alat, omega, ngm, ngl, igtongl, &
g, rhopol, nl, 1, gstart, gamma_only, vloc, &
force_environ )
!
IF ( env_static_permittivity .GT. 1.D0 ) THEN
!
IF ( TRIM(eps_mode) .NE. 'ionic' ) THEN
CALL force_lc( nat, tau, ityp, alat, omega, ngm, ngl, igtongl, &
g, rhopol, nl, 1, gstart, gamma_only, vloc, &
force_environ )
ELSE
CALL calc_fsolvent( dfftp%nnr, force_environ )
END IF
!
END IF
! ... Add the other environment contributions
!
CALL stop_clock('calc_fsolv')
CALL calc_fenviron( dfftp%nnr, nat, force_environ )
!
END IF
!
#endif
!
! Berry's phase electric field terms
!
......
......@@ -123,7 +123,8 @@ SUBROUTINE iosys()
tolrhopol_ => tolrhopol, &
env_surface_tension_ => env_surface_tension, &
delta_ => delta, &
env_pressure_ => env_pressure
env_pressure_ => env_pressure, &
env_slab_geometry, slab_axis
#endif
!
USE esm, ONLY: do_comp_esm, &
......@@ -1227,17 +1228,20 @@ SUBROUTINE iosys()
env_static_permittivity_ = 1.D0
env_surface_tension_ = 0.D0
env_pressure_ = 0.D0
env_slab_geometry = .false.
CASE ('water')
! water, experimental and SCCS tuned parameters
env_static_permittivity_ = 78.3D0
env_surface_tension_ = 50.D0*1.D-3*bohr_radius_si**2/rydberg_si
env_pressure_ = -0.35D0*1.D9/rydberg_si*bohr_radius_si**3
env_slab_geometry = .false.
CASE ('input')
! take values from input, this is the default option
env_static_permittivity_ = env_static_permittivity
env_surface_tension_ = &
env_surface_tension*1.D-3*bohr_radius_si**2/rydberg_si
env_pressure_ = env_pressure*1.D9/rydberg_si*bohr_radius_si**3
env_slab_geometry = .false.
CASE DEFAULT
call errore ('iosys','unrecognized value for environ_type',1)
END SELECT
......@@ -1289,6 +1293,35 @@ SUBROUTINE iosys()
do_comp_mt = .false.
do_makov_payne = .false.
!
#ifdef __ENVIRON
CASE( 'slabx' )
!
do_environ_ = .true.
env_slab_geometry = .true.
slab_axis = 1
do_makov_payne = .false.
do_comp_mt = .false.
do_comp_esm = .false.
!
CASE( 'slaby' )
!
do_environ_ = .true.
env_slab_geometry = .true.
slab_axis = 2
do_makov_payne = .false.
do_comp_mt = .false.
do_comp_esm = .false.
!
CASE( 'slabz' )
!
do_environ_ = .true.
env_slab_geometry = .true.
slab_axis = 3
do_makov_payne = .false.
do_comp_mt = .false.
do_comp_esm = .false.
!
#endif
CASE( 'none' )
!
do_makov_payne = .false.
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment