Add surface dipole correction [sghosh]

src/cp_control_types.F

@@ -468,6 +468,8 @@ MODULE cp_control_types
 
 ! **************************************************************************************************
 ! \brief Control parameters for a DFT calculation
+! \par History 
+!      10.2019 added variables related to surface dipole correction [Soumya Ghosh]
 ! **************************************************************************************************
    TYPE dft_control_type
       TYPE(admm_control_type), POINTER     :: admm_control
@@ -497,7 +499,8 @@ MODULE cp_control_types
                                               auto_basis_ri_xas = 1
       REAL(KIND=dp)                        :: relax_multiplicity, &
                                               sic_scaling_a, &
-                                              sic_scaling_b
+                                              sic_scaling_b, &
+                                              pos_dir_surf_dip
       LOGICAL                              :: do_tddfpt_calculation, &
                                               do_xas_calculation, &
                                               do_xas_tdp_calculation, &
@@ -523,6 +526,9 @@ MODULE cp_control_types
                                               apply_external_vxc, &
                                               read_external_vxc, &
                                               correct_surf_dip, &
+                                              surf_dip_correct_switch, &
+                                              switch_surf_dip, &
+                                              correct_el_density_dip, &
                                               do_sccs, &
                                               apply_embed_pot, &
                                               apply_dmfet_pot

src/cp_control_utils.F

@@ -322,6 +322,23 @@ SUBROUTINE read_dft_control(dft_control, dft_section)
       dft_control%correct_surf_dip = .FALSE.
       CALL section_vals_val_get(dft_section, "SURFACE_DIPOLE_CORRECTION", l_val=dft_control%correct_surf_dip)
       CALL section_vals_val_get(dft_section, "SURF_DIP_DIR", i_val=dft_control%dir_surf_dip)
+      dft_control%pos_dir_surf_dip = -1.0_dp
+      CALL section_vals_val_get(dft_section, "SURF_DIP_POS", r_val=dft_control%pos_dir_surf_dip)
+! another logical variable, surf_dip_correct_switch, is introduced for
+! implementation of "SURF_DIP_SWITCH" [SGh]
+      dft_control%switch_surf_dip = .FALSE.
+      dft_control%surf_dip_correct_switch = dft_control%correct_surf_dip
+      CALL section_vals_val_get(dft_section, "SURF_DIP_SWITCH", l_val=dft_control%switch_surf_dip)
+      dft_control%correct_el_density_dip = .FALSE.
+      CALL section_vals_val_get(dft_section, "CORE_CORR_DIP", l_val=dft_control%correct_el_density_dip)
+      IF (dft_control%correct_el_density_dip) THEN
+         IF (dft_control%correct_surf_dip) THEN
+            ! Do nothing, move on
+         ELSE
+            dft_control%correct_el_density_dip = .FALSE.
+            CPWARN("CORE_CORR_DIP keyword is activated only if SURFACE_DIPOLE_CORRECTION is TRUE")
+         ENDIF
+      ENDIF
 
       CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
                                 c_val=basis_set_file_name)

src/input_cp2k_dft.F

@@ -326,6 +326,41 @@ SUBROUTINE create_dft_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, &
+                          name="SURF_DIP_POS", &
+                          description="This keyword assigns an user defined position in Angstroms"// &
+                          " in the direction normal to the surface (given by SURF_DIP_DIR). "// &
+                          "The default value is -1.0_dp which appplies the correction at a position "// &
+                          " that has minimum electron density on the grid.", &
+                          usage="SURF_DIP_POS -1.0_dp", &
+                          default_r_val=-1.0_dp)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, &
+                          name="SURF_DIP_SWITCH", &
+                          description="WARNING: Experimental feature under development that will help the "// &
+                          "user to switch parameters to facilitate  SCF convergence. In its current form the "// &
+                          "surface dipole correction is switched off if the calculation does not converge in "// &
+                          "(0.5*MAX_SCF + 1) outer_scf steps. "// &
+                          "The default value is .FALSE.", &
+                          usage="SURF_DIP_SWITCH .TRUE.", &
+                          default_l_val=.FALSE., &
+                          lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, &
+                          name="CORE_CORR_DIP", &
+                          description="If the total CORE_CORRECTION is non-zero and surface dipole "// &
+                          "correction is switched on, presence of this keyword will adjust electron "// &
+                          "density via MO occupation to reflect the total CORE_CORRECTION. "// &
+                          "The default value is .FALSE.", &
+                          usage="CORE_CORR_DIP .TRUE.", &
+                          default_l_val=.FALSE., &
+                          lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
       NULLIFY (subsection)
       CALL create_scf_section(subsection)
       CALL section_add_subsection(section, subsection)

src/motion/pint_methods.F

@@ -2495,7 +2495,7 @@ SUBROUTINE pint_calc_energy(pint_env)
       CALL pint_calc_e_kin_beads_u(pint_env)
       CALL pint_calc_e_vir(pint_env)
 
-      CALL pint_calc_uf_h(pint_env, e_h)
+      CALL pint_calc_uf_h(pint_env, e_h=e_h)
       pint_env%e_pot_h = e_h
 
       SELECT CASE (pint_env%pimd_thermostat)

src/particle_methods.F

@@ -584,7 +584,7 @@ SUBROUTINE write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_secti
                                  z=z)
             CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
             WRITE (UNIT=iw, &
-                   FMT="(T2,I7,1X,I5,1X,A2,1X,I3,3F12.6,4X,F6.2,2X,F11.4)") &
+                   FMT="(T2,I7,1X,I5,1X,A2,1X,I3,3F12.6,4X,F8.4,2X,F11.4)") &
                iatom, ikind, element_symbol, z, &
                particle_set(iatom)%r(1:3)*conv, zeff, mass/massunit
          END DO

src/qs_environment.F

@@ -509,6 +509,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
                                                             lribas, orb_gradient, was_present
       REAL(dp)                                           :: alpha, ccore, ewald_rcut, maxocc, &
                                                             verlet_skin, zeff_correction
+      REAL(KIND=dp)                                      :: total_zeff_corr
       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
       TYPE(cp_logger_type), POINTER                      :: logger
       TYPE(dft_control_type), POINTER                    :: dft_control
@@ -523,7 +524,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
                                                             tmp_basis_set
       TYPE(local_rho_type), POINTER                      :: local_rho_set
       TYPE(lri_environment_type), POINTER                :: lri_env
-      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos, mos_aux_fit
+      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos, mos_aux_fit, mos_last_converged
       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
       TYPE(nddo_mpole_type), POINTER                     :: se_nddo_mpole
@@ -572,6 +573,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
       NULLIFY (dft_section)
       NULLIFY (et_coupling_section)
       NULLIFY (ks_env)
+      NULLIFY (mos_last_converged)
       dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
       qs_section => section_vals_get_subs_vals(dft_section, "QS")
       et_coupling_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%ET_COUPLING")
@@ -1298,6 +1300,11 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
 
       END IF
 
+      ! Read the total_zeff_corr here [SGh]
+      CALL get_qs_kind_set(qs_kind_set, total_zeff_corr=total_zeff_corr)
+      ! store it in qs_env
+      qs_env%total_zeff_corr = total_zeff_corr
+
       ! store the number of electrons once an for all
       CALL qs_subsys_set(subsys, &
                          nelectron_total=nelectron, &
@@ -1435,6 +1442,22 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
          CALL set_qs_env(qs_env, mos_aux_fit=mos_aux_fit)
       END IF
 
+! allocate mos when switch_surf_dip is triggered [SGh]
+      IF (dft_control%switch_surf_dip) THEN
+         ALLOCATE (mos_last_converged(dft_control%nspins))
+         DO ispin = 1, dft_control%nspins
+            NULLIFY (mos_last_converged(ispin)%mo_set)
+            CALL allocate_mo_set(mo_set=mos_last_converged(ispin)%mo_set, &
+                                 nao=n_ao, &
+                                 nmo=n_mo(ispin), &
+                                 nelectron=nelectron_spin(ispin), &
+                                 n_el_f=REAL(nelectron_spin(ispin), dp), &
+                                 maxocc=maxocc, &
+                                 flexible_electron_count=dft_control%relax_multiplicity)
+         END DO
+         CALL set_qs_env(qs_env, mos_last_converged=mos_last_converged)
+      ENDIF
+
       IF (.NOT. be_silent) THEN
          ! Print the DFT control parameters
          CALL write_dft_control(dft_control, dft_section)

src/qs_environment_types.F

@@ -312,6 +312,9 @@ MODULE qs_environment_types
       TYPE(polar_env_type), POINTER                         :: polar_env
       ! Resp charges
       REAL(KIND=dp), DIMENSION(:), POINTER                  :: rhs => NULL()
+      REAL(KIND=dp)                                         :: total_zeff_corr, surface_dipole_moment
+      LOGICAL                                               :: surface_dipole_switch_off
+      TYPE(mo_set_p_type), DIMENSION(:), POINTER            :: mos_last_converged
    END TYPE qs_environment_type
 
 ! **************************************************************************************************
@@ -492,6 +495,7 @@ MODULE qs_environment_types
 !> \param embed_pot ...
 !> \param spin_embed_pot ...
 !> \param polar_env ...
+!> \param mos_last_converged ... [SGh]
 !> \param rhs ...
 !> \date    23.01.2002
 !> \author  MK
@@ -523,7 +527,7 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce
                          lri_env, lri_density, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, &
                          mp2_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, &
                          s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, &
-                         gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, rhs)
+                         gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
       TYPE(qs_environment_type), POINTER                 :: qs_env
       TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, &
          POINTER                                         :: atomic_kind_set
@@ -657,6 +661,8 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce
       REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: gradient_history, variable_history
       TYPE(pw_p_type), OPTIONAL, POINTER                 :: embed_pot, spin_embed_pot
       TYPE(polar_env_type), OPTIONAL, POINTER            :: polar_env
+      TYPE(mo_set_p_type), DIMENSION(:), OPTIONAL, &
+         POINTER                                         :: mos_last_converged
       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: rhs
 
       CHARACTER(len=*), PARAMETER :: routineN = 'get_qs_env', routineP = moduleN//':'//routineN
@@ -679,6 +685,7 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce
       IF (PRESENT(qmmm_periodic)) qmmm_periodic = qs_env%qmmm_periodic
       IF (PRESENT(mos)) mos => qs_env%mos
       IF (PRESENT(mos_aux_fit)) mos_aux_fit => qs_env%mos_aux_fit
+      IF (PRESENT(mos_last_converged)) mos_last_converged => qs_env%mos_last_converged
       IF (PRESENT(ewald_env)) ewald_env => qs_env%ewald_env
       IF (PRESENT(ewald_pw)) ewald_pw => qs_env%ewald_pw
       IF (PRESENT(mpools)) mpools => qs_env%mpools
@@ -883,6 +890,7 @@ SUBROUTINE init_qs_env(qs_env, globenv)
       NULLIFY (qs_env%super_cell)
       NULLIFY (qs_env%mos)
       NULLIFY (qs_env%mos_aux_fit)
+      NULLIFY (qs_env%mos_last_converged)
       NULLIFY (qs_env%mpools)
       NULLIFY (qs_env%mpools_aux_fit)
       NULLIFY (qs_env%ewald_env)
@@ -967,6 +975,11 @@ SUBROUTINE init_qs_env(qs_env, globenv)
 
       qs_env%sim_time = 0._dp
       qs_env%sim_step = 0
+
+      qs_env%total_zeff_corr = 0.0_dp
+      qs_env%surface_dipole_moment = 0.0_dp
+      qs_env%surface_dipole_switch_off = .FALSE.
+
       ! Zero all variables containing results
       NULLIFY (qs_env%mo_derivs)
       NULLIFY (qs_env%mo_derivs_aux_fit)
@@ -1046,6 +1059,7 @@ END SUBROUTINE init_qs_env
 !> \param embed_pot ...
 !> \param spin_embed_pot ...
 !> \param polar_env ...
+!> \param mos_last_converged ... [SGh]
 !> \param rhs ...
 !> \date    23.01.2002
 !> \author  MK
@@ -1066,7 +1080,7 @@ SUBROUTINE set_qs_env(qs_env, super_cell, &
                          do_transport, transport_env, lri_env, lri_density, ec_env, dispersion_env, &
                          gcp_env, mp2_env, kg_env, force, &
                          kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, &
-                         spin_embed_pot, polar_env, rhs)
+                         spin_embed_pot, polar_env, mos_last_converged, rhs)
 
       TYPE(qs_environment_type), POINTER                 :: qs_env
       TYPE(cell_type), OPTIONAL, POINTER                 :: super_cell
@@ -1139,6 +1153,8 @@ SUBROUTINE set_qs_env(qs_env, super_cell, &
       REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: gradient_history, variable_history
       TYPE(pw_p_type), OPTIONAL, POINTER                 :: embed_pot, spin_embed_pot
       TYPE(polar_env_type), OPTIONAL, POINTER            :: polar_env
+      TYPE(mo_set_p_type), DIMENSION(:), OPTIONAL, &
+         POINTER                                         :: mos_last_converged
       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: rhs
 
       CHARACTER(len=*), PARAMETER :: routineN = 'set_qs_env', routineP = moduleN//':'//routineN
@@ -1160,6 +1176,7 @@ SUBROUTINE set_qs_env(qs_env, super_cell, &
       IF (PRESENT(qmmm_periodic)) qs_env%qmmm_periodic = qmmm_periodic
       IF (PRESENT(mos)) qs_env%mos => mos
       IF (PRESENT(mos_aux_fit)) qs_env%mos_aux_fit => mos_aux_fit
+      IF (PRESENT(mos_last_converged)) qs_env%mos_last_converged => mos_last_converged
       IF (PRESENT(ls_scf_env)) qs_env%ls_scf_env => ls_scf_env
       IF (PRESENT(almo_scf_env)) qs_env%almo_scf_env => almo_scf_env
       IF (PRESENT(do_transport)) qs_env%do_transport = do_transport
@@ -1395,6 +1412,12 @@ SUBROUTINE qs_env_release(qs_env)
                END DO
                DEALLOCATE (qs_env%mos_aux_fit)
             END IF
+            IF (ASSOCIATED(qs_env%mos_last_converged)) THEN
+               DO i = 1, SIZE(qs_env%mos_last_converged)
+                  CALL deallocate_mo_set(qs_env%mos_last_converged(i)%mo_set)
+               END DO
+               DEALLOCATE (qs_env%mos_last_converged)
+            END IF
 
             IF (ASSOCIATED(qs_env%mo_derivs)) THEN
                DO I = 1, SIZE(qs_env%mo_derivs)