1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> 09.2004 created [tlaino] 9!> \author Teodoro Laino 10! ************************************************************************************************** 11MODULE qmmm_elpot 12 USE cell_types, ONLY: cell_type 13 USE cp_log_handling, ONLY: cp_get_default_logger,& 14 cp_logger_get_default_io_unit,& 15 cp_logger_type 16 USE cp_output_handling, ONLY: cp_p_file,& 17 cp_print_key_finished_output,& 18 cp_print_key_should_output,& 19 cp_print_key_unit_nr 20 USE input_constants, ONLY: do_qmmm_coulomb,& 21 do_qmmm_gauss,& 22 do_qmmm_pcharge,& 23 do_qmmm_swave 24 USE input_section_types, ONLY: section_vals_type 25 USE kinds, ONLY: default_path_length,& 26 default_string_length,& 27 dp 28 USE mathconstants, ONLY: rootpi 29 USE memory_utilities, ONLY: reallocate 30 USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,& 31 qmmm_gaussian_type 32 USE qmmm_types_low, ONLY: qmmm_Pot_Type,& 33 qmmm_pot_p_type 34#include "./base/base_uses.f90" 35 36 IMPLICIT NONE 37 PRIVATE 38 39 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 40 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_elpot' 41 PUBLIC :: qmmm_potential_init 42 43CONTAINS 44 45! ************************************************************************************************** 46!> \brief Initialize the QMMM potential stored on vector, 47!> according the qmmm_coupl_type 48!> \param qmmm_coupl_type ... 49!> \param mm_el_pot_radius ... 50!> \param potentials ... 51!> \param pgfs ... 52!> \param mm_cell ... 53!> \param compatibility ... 54!> \param print_section ... 55!> \par History 56!> 09.2004 created [tlaino] 57!> \author Teodoro Laino 58! ************************************************************************************************** 59 SUBROUTINE qmmm_potential_init(qmmm_coupl_type, mm_el_pot_radius, potentials, & 60 pgfs, mm_cell, compatibility, print_section) 61 INTEGER, INTENT(IN) :: qmmm_coupl_type 62 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_el_pot_radius 63 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials 64 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs 65 TYPE(cell_type), POINTER :: mm_cell 66 LOGICAL, INTENT(IN) :: compatibility 67 TYPE(section_vals_type), POINTER :: print_section 68 69 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_potential_init', & 70 routineP = moduleN//':'//routineN 71 REAL(KIND=dp), PARAMETER :: dx = 0.05_dp 72 73 CHARACTER(LEN=default_path_length) :: myFormat 74 CHARACTER(LEN=default_string_length) :: rc_s 75 INTEGER :: I, ig, ig_start, J, K, myind, ndim, Np, & 76 output_unit, unit_nr 77 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 78 LOGICAL :: found 79 REAL(KIND=dp) :: A, G, rc, Rmax, Rmin, t, t1, t2, x 80 REAL(KIND=dp), DIMENSION(:), POINTER :: radius 81 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pot0_2 82 TYPE(cp_logger_type), POINTER :: logger 83 TYPE(qmmm_gaussian_type), POINTER :: pgf 84 85 logger => cp_get_default_logger() 86 output_unit = cp_logger_get_default_io_unit(logger) 87 Rmin = 0.0_dp 88 Rmax = SQRT(mm_cell%hmat(1, 1)**2 + & 89 mm_cell%hmat(2, 2)**2 + & 90 mm_cell%hmat(3, 3)**2) 91 np = CEILING(rmax/dx) + 1 92 ! 93 ! Preprocessing 94 ! 95 IF (SIZE(mm_el_pot_radius) /= 0) THEN 96 ALLOCATE (radius(1)) 97 radius(1) = mm_el_pot_radius(1) 98 ELSE 99 ALLOCATE (radius(0)) 100 END IF 101 Loop_on_all_values: DO I = 2, SIZE(mm_el_pot_radius) 102 Found = .FALSE. 103 Loop_on_found_values: DO J = 1, SIZE(radius) 104 IF (mm_el_pot_radius(i) .EQ. radius(j)) THEN 105 Found = .TRUE. 106 EXIT Loop_on_found_values 107 END IF 108 END DO Loop_on_found_values 109 IF (.NOT. Found) THEN 110 Ndim = SIZE(radius) 111 Ndim = Ndim + 1 112 CALL REALLOCATE(radius, 1, Ndim) 113 radius(Ndim) = mm_el_pot_radius(i) 114 END IF 115 END DO Loop_on_all_values 116 ! 117 CPASSERT(.NOT. ASSOCIATED(potentials)) 118 ALLOCATE (potentials(SIZE(radius))) 119 120 Potential_Type: DO K = 1, SIZE(radius) 121 122 rc = radius(K) 123 ALLOCATE (potentials(K)%Pot) 124 SELECT CASE (qmmm_coupl_type) 125 CASE (do_qmmm_coulomb) 126 NULLIFY (pot0_2) 127 CASE (do_qmmm_pcharge) 128 NULLIFY (pot0_2) 129 CASE (do_qmmm_gauss, do_qmmm_swave) 130 ALLOCATE (pot0_2(2, np)) 131 END SELECT 132 133 SELECT CASE (qmmm_coupl_type) 134 CASE (do_qmmm_coulomb, do_qmmm_pcharge) 135 ! Do Nothing 136 CASE (do_qmmm_gauss, do_qmmm_swave) 137 IF (qmmm_coupl_type == do_qmmm_gauss) THEN 138 ! Smooth Coulomb Potential :: Erf(x/rc)/x 139 pot0_2(1, 1) = 2.0_dp/(rootpi*rc) 140 pot0_2(2, 1) = 0.0_dp 141 x = 0.0_dp 142 DO i = 2, np 143 x = x + dx 144 pot0_2(1, i) = erf(x/rc)/x 145 t = 2._dp/(rootpi*x*rc)*EXP(-(x/rc)**2) 146 pot0_2(2, i) = (t - pot0_2(1, i)/x)*dx 147 END DO 148 ELSEIF (qmmm_coupl_type == do_qmmm_swave) THEN 149 ! S-wave expansion :: 1/x - exp(-2*x/rc) * ( 1/x - 1/rc ) 150 pot0_2(1, 1) = 1.0_dp/rc 151 pot0_2(2, 1) = 0.0_dp 152 x = 0.0_dp 153 DO i = 2, np 154 x = x + dx 155 t = EXP(-2.0_dp*x/rc)/rc 156 pot0_2(1, i) = (1.0_dp - t*(rc + x))/x 157 pot0_2(2, i) = ((t*(rc**2 + 2.0_dp*rc*x + 2.0_dp*x**2)/rc - 1.0_dp)/x**2)*dx 158 END DO 159 END IF 160 pgf => pgfs(K)%pgf 161 CPASSERT(pgf%Elp_Radius == rc) 162 ig_start = 1 163 IF (compatibility .AND. (qmmm_coupl_type == do_qmmm_gauss)) ig_start = 2 164 DO Ig = ig_start, pgf%number_of_gaussians 165 A = pgf%Ak(Ig) 166 G = pgf%Gk(Ig) 167 pot0_2(1, 1) = pot0_2(1, 1) - A 168 x = 0.0_dp 169 DO i = 2, np 170 x = x + dx 171 t = EXP(-(x/G)**2)*A 172 t1 = 1/G**2 173 t2 = t1*t 174 pot0_2(1, i) = pot0_2(1, i) - t 175 pot0_2(2, i) = pot0_2(2, i) + 2.0_dp*x*t2*dx 176 END DO 177 END DO 178 179 ! Print info on the unidimensional MM electrostatic potential 180 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, "MM_POTENTIAL") & 181 , cp_p_file)) THEN 182 WRITE (rc_s, '(F6.3)') rc 183 unit_nr = cp_print_key_unit_nr(logger, print_section, "MM_POTENTIAL", & 184 extension="_rc="//TRIM(ADJUSTL(rc_s))//".data") 185 IF (unit_nr > 0) THEN 186 WRITE (unit_nr, '(A)') "# MM ELECTROSTATIC POTENTIAL - UNIDIMENSIONAL - ATOMIC UNITS" 187 WRITE (unit_nr, '(A,I5)') "# MM ELECTROSTATIC POTENTIAL - Nr. of Gaussians:", pgf%number_of_gaussians 188 WRITE (unit_nr, '(A,T10,A,T30,A,T300,A)') "#", "Xval", "Gaussians", "LongRange" 189 myFormat = "T10,F15.9,T30," 190 DO Ig = 1, pgf%number_of_gaussians 191 myind = INDEX(myFormat, " ") 192 WRITE (myFormat(myind:), '(A6)') "F12.9," 193 END DO 194 myind = INDEX(myFormat, " ") - 1 195 myFormat = myFormat(1:myind)//"T300,F15.9" 196 myind = INDEX(myFormat, " ") - 1 197 x = 0.0_dp 198 DO i = 1, np 199 WRITE (unit_nr, '('//myFormat(1:myind)//')') & 200 x, (EXP(-(x/pgf%Gk(Ig))**2)*pgf%Ak(Ig), Ig=1, pgf%number_of_gaussians), pot0_2(1, i) 201 x = x + dx 202 END DO 203 END IF 204 CALL cp_print_key_finished_output(unit_nr, logger, print_section, & 205 "MM_POTENTIAL") 206 END IF 207 CASE DEFAULT 208 DEALLOCATE (potentials(K)%Pot) 209 IF (output_unit > 0) WRITE (output_unit, '(A)') " QMMM Potential - Spline Interpolation - not Initialized!" 210 CYCLE Potential_Type 211 END SELECT 212 NULLIFY (mm_atom_index) 213 ALLOCATE (mm_atom_index(1)) 214 ! Build mm_atom_index List 215 DO J = 1, SIZE(mm_el_pot_radius) 216 IF (rc .EQ. mm_el_pot_radius(J)) THEN 217 Ndim = SIZE(mm_atom_index) 218 mm_atom_index(Ndim) = J 219 CALL reallocate(mm_atom_index, 1, Ndim + 1) 220 ENDIF 221 END DO 222 CALL reallocate(mm_atom_index, 1, Ndim) 223 224 NULLIFY (potentials(K)%Pot%pot0_2) 225 CALL qmmm_pot_type_create(potentials(K)%Pot, pot0_2=pot0_2, & 226 Rmax=Rmax, Rmin=Rmin, dx=dx, Rc=rc, npts=np, & 227 mm_atom_index=mm_atom_index) 228 229 END DO Potential_Type 230 DEALLOCATE (radius) 231 END SUBROUTINE qmmm_potential_init 232 233! ************************************************************************************************** 234!> \brief Creates the qmmm_pot_type structure 235!> \param Pot ... 236!> \param pot0_2 ... 237!> \param Rmax ... 238!> \param Rmin ... 239!> \param dx ... 240!> \param npts ... 241!> \param rc ... 242!> \param mm_atom_index ... 243!> \par History 244!> 09.2004 created [tlaino] 245!> \author Teodoro Laino 246! ************************************************************************************************** 247 SUBROUTINE qmmm_pot_type_create(Pot, pot0_2, Rmax, Rmin, dx, npts, rc, & 248 mm_atom_index) 249 TYPE(qmmm_Pot_Type), POINTER :: Pot 250 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pot0_2 251 REAL(KIND=dp), INTENT(IN) :: Rmax, Rmin, dx 252 INTEGER, INTENT(IN) :: npts 253 REAL(KIND=dp), INTENT(IN) :: Rc 254 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 255 256 Pot%pot0_2 => pot0_2 257 Pot%mm_atom_index => mm_atom_index 258 Pot%Rmax = Rmax 259 Pot%Rmin = Rmin 260 Pot%Rc = rc 261 Pot%dx = dx 262 Pot%npts = npts 263 264 END SUBROUTINE qmmm_pot_type_create 265 266END MODULE qmmm_elpot 267