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!> Subroutine input_torsions changed (DG) 05-Dec-2000 9!> Output formats changed (DG) 05-Dec-2000 10!> JGH (26-01-2002) : force field parameters stored in tables, not in 11!> matrices. Input changed to have parameters labeled by the position 12!> and not atom pairs (triples etc) 13!> Teo (11.2005) : Moved all information on force field pair_potential to 14!> a much lighter memory structure 15!> Teo 09.2006 : Split all routines force_field I/O in a separate file 16!> \author CJM 17! ************************************************************************************************** 18MODULE force_fields_input 19 USE bibliography, ONLY: Siepmann1995,& 20 Tersoff1988,& 21 Tosi1964a,& 22 Tosi1964b,& 23 Yamada2000,& 24 cite_reference 25 USE cp_log_handling, ONLY: cp_get_default_logger,& 26 cp_logger_type,& 27 cp_to_string 28 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 29 cp_print_key_unit_nr 30 USE cp_para_types, ONLY: cp_para_env_type 31 USE cp_parser_methods, ONLY: parser_get_next_line 32 USE cp_parser_types, ONLY: cp_parser_type,& 33 parser_create,& 34 parser_release 35 USE cp_units, ONLY: cp_unit_to_cp2k 36 USE damping_dipole_types, ONLY: damping_info_type 37 USE force_field_kind_types, ONLY: do_ff_amber,& 38 do_ff_charmm,& 39 do_ff_g87,& 40 do_ff_g96,& 41 do_ff_opls,& 42 do_ff_undef,& 43 legendre_data_type 44 USE force_field_types, ONLY: force_field_type,& 45 input_info_type 46 USE force_fields_util, ONLY: get_generic_info 47 USE input_section_types, ONLY: section_vals_get,& 48 section_vals_get_subs_vals,& 49 section_vals_type,& 50 section_vals_val_get 51 USE kinds, ONLY: default_string_length,& 52 dp 53 USE mathconstants, ONLY: pi 54 USE mathlib, ONLY: invert_matrix 55 USE memory_utilities, ONLY: reallocate 56 USE pair_potential_types, ONLY: & 57 b4_type, bm_type, do_potential_single_allocation, ea_type, eam_pot_type, ft_pot_type, & 58 ft_type, ftd_type, gp_type, gw_type, ip_type, ipbv_pot_type, lj_charmm_type, & 59 no_potential_single_allocation, pair_potential_p_type, pair_potential_reallocate, & 60 potential_single_allocation, quip_type, siepmann_type, tersoff_type, wl_type 61 USE shell_potential_types, ONLY: shell_p_create,& 62 shell_p_type 63 USE string_utilities, ONLY: uppercase 64#include "./base/base_uses.f90" 65 66 IMPLICIT NONE 67 68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_input' 69 70 PRIVATE 71 PUBLIC :: read_force_field_section, & 72 read_lj_section, & 73 read_wl_section, & 74 read_gd_section, & 75 read_gp_section, & 76 read_chrg_section 77 78CONTAINS 79 80! ************************************************************************************************** 81!> \brief Reads the force_field input section 82!> \param ff_section ... 83!> \param mm_section ... 84!> \param ff_type ... 85!> \param para_env ... 86!> \author teo 87! ************************************************************************************************** 88 SUBROUTINE read_force_field_section1(ff_section, mm_section, ff_type, para_env) 89 TYPE(section_vals_type), POINTER :: ff_section, mm_section 90 TYPE(force_field_type), INTENT(INOUT) :: ff_type 91 TYPE(cp_para_env_type), POINTER :: para_env 92 93 CHARACTER(len=*), PARAMETER :: routineN = 'read_force_field_section1', & 94 routineP = moduleN//':'//routineN 95 96 INTEGER :: nb4, nbends, nbm, nbmhft, nbmhftd, nbonds, nchg, neam, ngd, ngp, nimpr, nipbv, & 97 nlj, nopbend, nquip, nshell, nsiepmann, ntersoff, ntors, ntot, nubs, nwl 98 LOGICAL :: explicit, unique_spline 99 REAL(KIND=dp) :: min_eps_spline_allowed 100 TYPE(input_info_type), POINTER :: inp_info 101 TYPE(section_vals_type), POINTER :: tmp_section, tmp_section2 102 103 INTEGER::i 104 105 NULLIFY (tmp_section, tmp_section2) 106 inp_info => ff_type%inp_info 107 CALL section_vals_val_get(ff_section, "PARMTYPE", i_val=ff_type%ff_type) 108 CALL section_vals_val_get(ff_section, "EI_SCALE14", r_val=ff_type%ei_scale14) 109 CALL section_vals_val_get(ff_section, "VDW_SCALE14", r_val=ff_type%vdw_scale14) 110 CALL section_vals_val_get(ff_section, "SPLINE%RCUT_NB", r_val=ff_type%rcut_nb) 111 CALL section_vals_val_get(ff_section, "SPLINE%R0_NB", r_val=ff_type%rlow_nb) 112 CALL section_vals_val_get(ff_section, "SPLINE%EPS_SPLINE", r_val=ff_type%eps_spline) 113 CALL section_vals_val_get(ff_section, "SPLINE%EMAX_SPLINE", r_val=ff_type%emax_spline) 114 CALL section_vals_val_get(ff_section, "SPLINE%EMAX_ACCURACY", r_val=ff_type%max_energy) 115 CALL section_vals_val_get(ff_section, "SPLINE%NPOINTS", i_val=ff_type%npoints) 116 CALL section_vals_val_get(ff_section, "IGNORE_MISSING_CRITICAL_PARAMS", l_val=ff_type%ignore_missing_critical) 117 CPASSERT(ff_type%max_energy <= ff_type%emax_spline) 118 ! Read the parameter file name only if the force field type requires it.. 119 SELECT CASE (ff_type%ff_type) 120 CASE (do_ff_charmm, do_ff_amber, do_ff_g96, do_ff_g87) 121 CALL section_vals_val_get(ff_section, "PARM_FILE_NAME", c_val=ff_type%ff_file_name) 122 IF (TRIM(ff_type%ff_file_name) == "") & 123 CPABORT("Force Field Parameter's filename is empty! Please check your input file.") 124 CASE (do_ff_undef) 125 ! Do Nothing 126 CASE DEFAULT 127 CPABORT("Force field type not implemented") 128 END SELECT 129 ! Numerical Accuracy: 130 ! the factors here should depend on the energy and on the shape of each potential mapped 131 ! with splines. this would make everything un-necessarily complicated. Let's just be safe 132 ! and assume that we cannot achieve an accuracy on the spline 2 orders of magnitude more 133 ! than the smallest representable number (taking into account also the max_energy for the 134 ! spline generation 135 min_eps_spline_allowed = 20.0_dp*MAX(ff_type%max_energy, 10.0_dp)*EPSILON(0.0_dp) 136 IF (ff_type%eps_spline < min_eps_spline_allowed) THEN 137 CALL cp_warn(__LOCATION__, & 138 "Requested spline accuracy ("//TRIM(cp_to_string(ff_type%eps_spline))//" ) "// & 139 "is smaller than the minimum value allowed ("//TRIM(cp_to_string(min_eps_spline_allowed))// & 140 " ) with the present machine precision ("//TRIM(cp_to_string(EPSILON(0.0_dp)))//" ). "// & 141 "New EPS_SPLINE value ("//TRIM(cp_to_string(min_eps_spline_allowed))//" ). ") 142 ff_type%eps_spline = min_eps_spline_allowed 143 END IF 144 CALL section_vals_val_get(ff_section, "SHIFT_CUTOFF", l_val=ff_type%shift_cutoff) 145 CALL section_vals_val_get(ff_section, "SPLINE%UNIQUE_SPLINE", l_val=unique_spline) 146 ! Single spline 147 potential_single_allocation = no_potential_single_allocation 148 IF (unique_spline) potential_single_allocation = do_potential_single_allocation 149 150 CALL section_vals_val_get(ff_section, "MULTIPLE_POTENTIAL", l_val=ff_type%multiple_potential) 151 CALL section_vals_val_get(ff_section, "DO_NONBONDED", l_val=ff_type%do_nonbonded) 152 tmp_section => section_vals_get_subs_vals(ff_section, "NONBONDED") 153 CALL section_vals_get(tmp_section, explicit=explicit) 154 IF (explicit .AND. ff_type%do_nonbonded) THEN 155 tmp_section2 => section_vals_get_subs_vals(tmp_section, "LENNARD-JONES") 156 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nlj) 157 ntot = 0 158 IF (explicit) THEN 159 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nlj, lj_charmm=.TRUE.) 160 CALL read_lj_section(inp_info%nonbonded, tmp_section2, ntot) 161 END IF 162 163 tmp_section2 => section_vals_get_subs_vals(tmp_section, "WILLIAMS") 164 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nwl) 165 ntot = nlj 166 IF (explicit) THEN 167 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nwl, williams=.TRUE.) 168 CALL read_wl_section(inp_info%nonbonded, tmp_section2, ntot) 169 END IF 170 171 tmp_section2 => section_vals_get_subs_vals(tmp_section, "EAM") 172 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=neam) 173 ntot = nlj + nwl 174 IF (explicit) THEN 175 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + neam, eam=.TRUE.) 176 CALL read_eam_section(inp_info%nonbonded, tmp_section2, ntot, para_env, mm_section) 177 END IF 178 179 tmp_section2 => section_vals_get_subs_vals(tmp_section, "GOODWIN") 180 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngd) 181 ntot = nlj + nwl + neam 182 IF (explicit) THEN 183 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ngd, goodwin=.TRUE.) 184 CALL read_gd_section(inp_info%nonbonded, tmp_section2, ntot) 185 END IF 186 187 tmp_section2 => section_vals_get_subs_vals(tmp_section, "IPBV") 188 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nipbv) 189 ntot = nlj + nwl + neam + ngd 190 IF (explicit) THEN 191 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nipbv, ipbv=.TRUE.) 192 CALL read_ipbv_section(inp_info%nonbonded, tmp_section2, ntot) 193 END IF 194 195 tmp_section2 => section_vals_get_subs_vals(tmp_section, "BMHFT") 196 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nbmhft) 197 ntot = nlj + nwl + neam + ngd + nipbv 198 IF (explicit) THEN 199 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nbmhft, bmhft=.TRUE.) 200 CALL read_bmhft_section(inp_info%nonbonded, tmp_section2, ntot) 201 END IF 202 203 tmp_section2 => section_vals_get_subs_vals(tmp_section, "BMHFTD") 204 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nbmhftd) 205 ntot = nlj + nwl + neam + ngd + nipbv + nbmhft 206 IF (explicit) THEN 207 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nbmhftd, bmhftd=.TRUE.) 208 CALL read_bmhftd_section(inp_info%nonbonded, tmp_section2, ntot) 209 END IF 210 211 tmp_section2 => section_vals_get_subs_vals(tmp_section, "BUCK4RANGES") 212 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nb4) 213 ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd 214 IF (explicit) THEN 215 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nb4, buck4r=.TRUE.) 216 CALL read_b4_section(inp_info%nonbonded, tmp_section2, ntot) 217 END IF 218 219 tmp_section2 => section_vals_get_subs_vals(tmp_section, "BUCKMORSE") 220 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nbm) 221 ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 222 IF (explicit) THEN 223 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nbm, buckmo=.TRUE.) 224 CALL read_bm_section(inp_info%nonbonded, tmp_section2, ntot) 225 END IF 226 227 tmp_section2 => section_vals_get_subs_vals(tmp_section, "GENPOT") 228 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngp) 229 ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm 230 IF (explicit) THEN 231 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ngp, gp=.TRUE.) 232 CALL read_gp_section(inp_info%nonbonded, tmp_section2, ntot) 233 END IF 234 tmp_section2 => section_vals_get_subs_vals(tmp_section, "TERSOFF") 235 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ntersoff) 236 ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp 237 IF (explicit) THEN 238 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ntersoff, tersoff=.TRUE.) 239 CALL read_tersoff_section(inp_info%nonbonded, tmp_section2, ntot, tmp_section2) 240 END IF 241 tmp_section2 => section_vals_get_subs_vals(tmp_section, "SIEPMANN") 242 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nsiepmann) 243 ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff 244 IF (explicit) THEN 245 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nsiepmann, siepmann=.TRUE.) 246 CALL read_siepmann_section(inp_info%nonbonded, tmp_section2, ntot, tmp_section2) 247 END IF 248 249 tmp_section2 => section_vals_get_subs_vals(tmp_section, "quip") 250 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nquip) 251 ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + nsiepmann 252 IF (explicit) THEN 253 CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nquip, quip=.TRUE.) 254 CALL read_quip_section(inp_info%nonbonded, tmp_section2, ntot) 255 END IF 256 END IF 257 258 tmp_section => section_vals_get_subs_vals(ff_section, "NONBONDED14") 259 CALL section_vals_get(tmp_section, explicit=explicit) 260 IF (explicit .AND. ff_type%do_nonbonded) THEN 261 tmp_section2 => section_vals_get_subs_vals(tmp_section, "LENNARD-JONES") 262 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nlj) 263 ntot = 0 264 IF (explicit) THEN 265 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, ntot + nlj, lj_charmm=.TRUE.) 266 CALL read_lj_section(inp_info%nonbonded14, tmp_section2, ntot) 267 END IF 268 tmp_section2 => section_vals_get_subs_vals(tmp_section, "WILLIAMS") 269 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nwl) 270 ntot = nlj 271 IF (explicit) THEN 272 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, ntot + nwl, williams=.TRUE.) 273 CALL read_wl_section(inp_info%nonbonded14, tmp_section2, ntot) 274 END IF 275 tmp_section2 => section_vals_get_subs_vals(tmp_section, "GOODWIN") 276 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngd) 277 ntot = nlj + nwl 278 IF (explicit) THEN 279 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, ntot + ngd, goodwin=.TRUE.) 280 CALL read_gd_section(inp_info%nonbonded14, tmp_section2, ntot) 281 END IF 282 tmp_section2 => section_vals_get_subs_vals(tmp_section, "GENPOT") 283 CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngp) 284 ntot = nlj + nwl + ngd 285 IF (explicit) THEN 286 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, ntot + ngp, gp=.TRUE.) 287 CALL read_gp_section(inp_info%nonbonded14, tmp_section2, ntot) 288 END IF 289 END IF 290 291 tmp_section => section_vals_get_subs_vals(ff_section, "CHARGE") 292 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nchg) 293 IF (explicit) THEN 294 ntot = 0 295 CALL reallocate(inp_info%charge_atm, 1, nchg) 296 CALL reallocate(inp_info%charge, 1, nchg) 297 CALL read_chrg_section(inp_info%charge_atm, inp_info%charge, tmp_section, ntot) 298 END IF 299 tmp_section => section_vals_get_subs_vals(ff_section, "DIPOLE") 300 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nchg) 301 IF (explicit) THEN 302 ntot = 0 303 CALL reallocate(inp_info%apol_atm, 1, nchg) 304 CALL reallocate(inp_info%apol, 1, nchg) 305 CALL read_apol_section(inp_info%apol_atm, inp_info%apol, inp_info%damping_list, & 306 tmp_section, ntot) 307 END IF 308 tmp_section => section_vals_get_subs_vals(ff_section, "QUADRUPOLE") 309 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nchg) 310 IF (explicit) THEN 311 ntot = 0 312 CALL reallocate(inp_info%cpol_atm, 1, nchg) 313 CALL reallocate(inp_info%cpol, 1, nchg) 314 CALL read_cpol_section(inp_info%cpol_atm, inp_info%cpol, tmp_section, ntot) 315 END IF 316 tmp_section => section_vals_get_subs_vals(ff_section, "SHELL") 317 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nshell) 318 IF (explicit) THEN 319 ntot = 0 320 CALL shell_p_create(inp_info%shell_list, nshell) 321 CALL read_shell_section(inp_info%shell_list, tmp_section, ntot) 322 END IF 323 324 tmp_section => section_vals_get_subs_vals(ff_section, "BOND") 325 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nbonds) 326 IF (explicit) THEN 327 ntot = 0 328 CALL reallocate(inp_info%bond_kind, 1, nbonds) 329 CALL reallocate(inp_info%bond_a, 1, nbonds) 330 CALL reallocate(inp_info%bond_b, 1, nbonds) 331 CALL reallocate(inp_info%bond_k, 1, 3, 1, nbonds) 332 CALL reallocate(inp_info%bond_r0, 1, nbonds) 333 CALL reallocate(inp_info%bond_cs, 1, nbonds) 334 CALL read_bonds_section(inp_info%bond_kind, inp_info%bond_a, inp_info%bond_b, inp_info%bond_k, & 335 inp_info%bond_r0, inp_info%bond_cs, tmp_section, ntot) 336 END IF 337 tmp_section => section_vals_get_subs_vals(ff_section, "BEND") 338 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nbends) 339 IF (explicit) THEN 340 ntot = 0 341 CALL reallocate(inp_info%bend_kind, 1, nbends) 342 CALL reallocate(inp_info%bend_a, 1, nbends) 343 CALL reallocate(inp_info%bend_b, 1, nbends) 344 CALL reallocate(inp_info%bend_c, 1, nbends) 345 CALL reallocate(inp_info%bend_k, 1, nbends) 346 CALL reallocate(inp_info%bend_theta0, 1, nbends) 347 CALL reallocate(inp_info%bend_cb, 1, nbends) 348 CALL reallocate(inp_info%bend_r012, 1, nbends) 349 CALL reallocate(inp_info%bend_r032, 1, nbends) 350 CALL reallocate(inp_info%bend_kbs12, 1, nbends) 351 CALL reallocate(inp_info%bend_kbs32, 1, nbends) 352 CALL reallocate(inp_info%bend_kss, 1, nbends) 353 IF (ASSOCIATED(inp_info%bend_legendre)) THEN 354 DO i = 1, SIZE(inp_info%bend_legendre) 355 IF (ASSOCIATED(inp_info%bend_legendre(i)%coeffs)) THEN 356 DEALLOCATE (inp_info%bend_legendre(i)%coeffs) 357 NULLIFY (inp_info%bend_legendre(i)%coeffs) 358 END IF 359 END DO 360 DEALLOCATE (inp_info%bend_legendre) 361 NULLIFY (inp_info%bend_legendre) 362 END IF 363 ALLOCATE (inp_info%bend_legendre(1:nbends)) 364 DO i = 1, SIZE(inp_info%bend_legendre(1:nbends)) 365 NULLIFY (inp_info%bend_legendre(i)%coeffs) 366 inp_info%bend_legendre(i)%order = 0 367 END DO 368 CALL read_bends_section(inp_info%bend_kind, inp_info%bend_a, inp_info%bend_b, inp_info%bend_c, & 369 inp_info%bend_k, inp_info%bend_theta0, inp_info%bend_cb, & 370 inp_info%bend_r012, inp_info%bend_r032, inp_info%bend_kbs12, & 371 inp_info%bend_kbs32, inp_info%bend_kss, & 372 inp_info%bend_legendre, tmp_section, ntot) 373 END IF 374 tmp_section => section_vals_get_subs_vals(ff_section, "BEND") 375 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nubs) 376 IF (explicit) THEN 377 ntot = 0 378 CALL reallocate(inp_info%ub_kind, 1, nubs) 379 CALL reallocate(inp_info%ub_a, 1, nubs) 380 CALL reallocate(inp_info%ub_b, 1, nubs) 381 CALL reallocate(inp_info%ub_c, 1, nubs) 382 CALL reallocate(inp_info%ub_k, 1, 3, 1, nubs) 383 CALL reallocate(inp_info%ub_r0, 1, nubs) 384 CALL read_ubs_section(inp_info%ub_kind, inp_info%ub_a, inp_info%ub_b, inp_info%ub_c, & 385 inp_info%ub_k, inp_info%ub_r0, tmp_section, ntot) 386 END IF 387 tmp_section => section_vals_get_subs_vals(ff_section, "TORSION") 388 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=ntors) 389 IF (explicit) THEN 390 ntot = 0 391 CALL reallocate(inp_info%torsion_kind, 1, ntors) 392 CALL reallocate(inp_info%torsion_a, 1, ntors) 393 CALL reallocate(inp_info%torsion_b, 1, ntors) 394 CALL reallocate(inp_info%torsion_c, 1, ntors) 395 CALL reallocate(inp_info%torsion_d, 1, ntors) 396 CALL reallocate(inp_info%torsion_k, 1, ntors) 397 CALL reallocate(inp_info%torsion_m, 1, ntors) 398 CALL reallocate(inp_info%torsion_phi0, 1, ntors) 399 CALL read_torsions_section(inp_info%torsion_kind, inp_info%torsion_a, inp_info%torsion_b, & 400 inp_info%torsion_c, inp_info%torsion_d, inp_info%torsion_k, inp_info%torsion_phi0, & 401 inp_info%torsion_m, tmp_section, ntot) 402 END IF 403 404 tmp_section => section_vals_get_subs_vals(ff_section, "IMPROPER") 405 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nimpr) 406 IF (explicit) THEN 407 ntot = 0 408 CALL reallocate(inp_info%impr_kind, 1, nimpr) 409 CALL reallocate(inp_info%impr_a, 1, nimpr) 410 CALL reallocate(inp_info%impr_b, 1, nimpr) 411 CALL reallocate(inp_info%impr_c, 1, nimpr) 412 CALL reallocate(inp_info%impr_d, 1, nimpr) 413 CALL reallocate(inp_info%impr_k, 1, nimpr) 414 CALL reallocate(inp_info%impr_phi0, 1, nimpr) 415 CALL read_improper_section(inp_info%impr_kind, inp_info%impr_a, inp_info%impr_b, & 416 inp_info%impr_c, inp_info%impr_d, inp_info%impr_k, inp_info%impr_phi0, & 417 tmp_section, ntot) 418 END IF 419 420 tmp_section => section_vals_get_subs_vals(ff_section, "OPBEND") 421 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nopbend) 422 IF (explicit) THEN 423 ntot = 0 424 CALL reallocate(inp_info%opbend_kind, 1, nopbend) 425 CALL reallocate(inp_info%opbend_a, 1, nopbend) 426 CALL reallocate(inp_info%opbend_b, 1, nopbend) 427 CALL reallocate(inp_info%opbend_c, 1, nopbend) 428 CALL reallocate(inp_info%opbend_d, 1, nopbend) 429 CALL reallocate(inp_info%opbend_k, 1, nopbend) 430 CALL reallocate(inp_info%opbend_phi0, 1, nopbend) 431 CALL read_opbend_section(inp_info%opbend_kind, inp_info%opbend_a, inp_info%opbend_b, & 432 inp_info%opbend_c, inp_info%opbend_d, inp_info%opbend_k, inp_info%opbend_phi0, & 433 tmp_section, ntot) 434 END IF 435 436 END SUBROUTINE read_force_field_section1 437 438! ************************************************************************************************** 439!> \brief Set up of the IPBV force fields 440!> \param at1 ... 441!> \param at2 ... 442!> \param ipbv ... 443!> \author teo 444! ************************************************************************************************** 445 SUBROUTINE set_IPBV_ff(at1, at2, ipbv) 446 CHARACTER(LEN=*), INTENT(IN) :: at1, at2 447 TYPE(ipbv_pot_type), POINTER :: ipbv 448 449 CHARACTER(len=*), PARAMETER :: routineN = 'set_IPBV_ff', routineP = moduleN//':'//routineN 450 451 IF ((at1(1:1) == 'O') .AND. (at2(1:1) == 'O')) THEN 452 ipbv%rcore = 0.9_dp ! a.u. 453 ipbv%m = -1.2226442563398141E+11_dp ! Kelvin/a.u. 454 ipbv%b = 1.1791292385486696E+11_dp ! Hartree 455 456 ! Hartree*a.u.^2 457 ipbv%a(2) = 4.786380682394_dp 458 ipbv%a(3) = -1543.407053545_dp 459 ipbv%a(4) = 88783.31188529_dp 460 ipbv%a(5) = -2361200.155376_dp 461 ipbv%a(6) = 35940504.84679_dp 462 ipbv%a(7) = -339762743.6358_dp 463 ipbv%a(8) = 2043874926.466_dp 464 ipbv%a(9) = -7654856796.383_dp 465 ipbv%a(10) = 16195251405.65_dp 466 ipbv%a(11) = -13140392992.18_dp 467 ipbv%a(12) = -9285572894.245_dp 468 ipbv%a(13) = 8756947519.029_dp 469 ipbv%a(14) = 15793297761.67_dp 470 ipbv%a(15) = 12917180227.21_dp 471 ELSEIF (((at1(1:1) == 'O') .AND. (at2(1:1) == 'H')) .OR. & 472 ((at1(1:1) == 'H') .AND. (at2(1:1) == 'O'))) THEN 473 ipbv%rcore = 2.95_dp ! a.u. 474 475 ipbv%m = -0.004025691139759147_dp ! Hartree/a.u. 476 ipbv%b = -2.193731138097428_dp ! Hartree 477 ! Hartree*a.u.^2 478 ipbv%a(2) = -195.7716013277_dp 479 ipbv%a(3) = 15343.78613395_dp 480 ipbv%a(4) = -530864.4586516_dp 481 ipbv%a(5) = 10707934.39058_dp 482 ipbv%a(6) = -140099704.7890_dp 483 ipbv%a(7) = 1250943273.785_dp 484 ipbv%a(8) = -7795458330.676_dp 485 ipbv%a(9) = 33955897217.31_dp 486 ipbv%a(10) = -101135640744.0_dp 487 ipbv%a(11) = 193107995718.7_dp 488 ipbv%a(12) = -193440560940.0_dp 489 ipbv%a(13) = -4224406093.918E0_dp 490 ipbv%a(14) = 217192386506.5E0_dp 491 ipbv%a(15) = -157581228915.5_dp 492 ELSEIF ((at1(1:1) == 'H') .AND. (at2(1:1) == 'H')) THEN 493 ipbv%rcore = 3.165_dp ! a.u. 494 ipbv%m = 0.002639704108787555_dp ! Hartree/a.u. 495 ipbv%b = -0.2735482611857583_dp ! Hartree 496 ! Hartree*a.u.^2 497 ipbv%a(2) = -26.29456010782_dp 498 ipbv%a(3) = 2373.352548248_dp 499 ipbv%a(4) = -93880.43551360_dp 500 ipbv%a(5) = 2154624.884809_dp 501 ipbv%a(6) = -31965151.34955_dp 502 ipbv%a(7) = 322781785.3278_dp 503 ipbv%a(8) = -2271097368.668_dp 504 ipbv%a(9) = 11169163192.90_dp 505 ipbv%a(10) = -37684457778.47_dp 506 ipbv%a(11) = 82562104256.03_dp 507 ipbv%a(12) = -100510435213.4_dp 508 ipbv%a(13) = 24570342714.65E0_dp 509 ipbv%a(14) = 88766181532.94E0_dp 510 ipbv%a(15) = -79705131323.98_dp 511 ELSE 512 CPABORT("IPBV only for WATER") 513 ENDIF 514 END SUBROUTINE set_IPBV_ff 515 516! ************************************************************************************************** 517!> \brief Set up of the BMHFT force fields 518!> \param at1 ... 519!> \param at2 ... 520!> \param ft ... 521!> \author teo 522! ************************************************************************************************** 523 SUBROUTINE set_BMHFT_ff(at1, at2, ft) 524 CHARACTER(LEN=*), INTENT(IN) :: at1, at2 525 TYPE(ft_pot_type), POINTER :: ft 526 527 CHARACTER(len=*), PARAMETER :: routineN = 'set_BMHFT_ff', routineP = moduleN//':'//routineN 528 529 ft%b = cp_unit_to_cp2k(3.1545_dp, "angstrom^-1") 530 IF ((at1(1:2) == 'NA') .AND. (at2(1:2) == 'NA')) THEN 531 ft%a = cp_unit_to_cp2k(424.097_dp, "eV") 532 ft%c = cp_unit_to_cp2k(1.05_dp, "eV*angstrom^6") 533 ft%d = cp_unit_to_cp2k(0.499_dp, "eV*angstrom^8") 534 ELSEIF (((at1(1:2) == 'NA') .AND. (at2(1:2) == 'CL')) .OR. & 535 ((at1(1:2) == 'CL') .AND. (at2(1:2) == 'NA'))) THEN 536 ft%a = cp_unit_to_cp2k(1256.31_dp, "eV") 537 ft%c = cp_unit_to_cp2k(7.00_dp, "eV*angstrom^6") 538 ft%d = cp_unit_to_cp2k(8.676_dp, "eV*angstrom^8") 539 ELSEIF ((at1(1:2) == 'CL') .AND. (at2(1:2) == 'CL')) THEN 540 ft%a = cp_unit_to_cp2k(3488.998_dp, "eV") 541 ft%c = cp_unit_to_cp2k(72.50_dp, "eV*angstrom^6") 542 ft%d = cp_unit_to_cp2k(145.427_dp, "eV*angstrom^8") 543 ELSE 544 CPABORT("BMHFT only for NaCl") 545 ENDIF 546 547 END SUBROUTINE set_BMHFT_ff 548 549! ************************************************************************************************** 550!> \brief Set up of the BMHFTD force fields 551!> \author Mathieu Salanne 05.2010 552! ************************************************************************************************** 553 SUBROUTINE set_BMHFTD_ff() 554 555 CHARACTER(len=*), PARAMETER :: routineN = 'set_BMHFTD_ff', routineP = moduleN//':'//routineN 556 557 CPABORT("No default parameters present for BMHFTD") 558 559 END SUBROUTINE set_BMHFTD_ff 560 561! ************************************************************************************************** 562!> \brief Reads the EAM section 563!> \param nonbonded ... 564!> \param section ... 565!> \param start ... 566!> \param para_env ... 567!> \param mm_section ... 568!> \author teo 569! ************************************************************************************************** 570 SUBROUTINE read_eam_section(nonbonded, section, start, para_env, mm_section) 571 TYPE(pair_potential_p_type), POINTER :: nonbonded 572 TYPE(section_vals_type), POINTER :: section 573 INTEGER, INTENT(IN) :: start 574 TYPE(cp_para_env_type), POINTER :: para_env 575 TYPE(section_vals_type), POINTER :: mm_section 576 577 CHARACTER(len=*), PARAMETER :: routineN = 'read_eam_section', & 578 routineP = moduleN//':'//routineN 579 580 CHARACTER(LEN=default_string_length), & 581 DIMENSION(:), POINTER :: atm_names 582 INTEGER :: isec, n_items 583 584 CALL section_vals_get(section, n_repetition=n_items) 585 DO isec = 1, n_items 586 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 587 588 nonbonded%pot(start + isec)%pot%type = ea_type 589 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 590 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 591 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 592 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 593 CALL section_vals_val_get(section, "PARM_FILE_NAME", i_rep_section=isec, & 594 c_val=nonbonded%pot(start + isec)%pot%set(1)%eam%eam_file_name) 595 CALL read_eam_data(nonbonded%pot(start + isec)%pot%set(1)%eam, para_env, mm_section) 596 nonbonded%pot(start + isec)%pot%rcutsq = nonbonded%pot(start + isec)%pot%set(1)%eam%acutal**2 597 END DO 598 END SUBROUTINE read_eam_section 599 600! ************************************************************************************************** 601!> \brief Reads the QUIP section 602!> \param nonbonded ... 603!> \param section ... 604!> \param start ... 605!> \author teo 606! ************************************************************************************************** 607 SUBROUTINE read_quip_section(nonbonded, section, start) 608 TYPE(pair_potential_p_type), POINTER :: nonbonded 609 TYPE(section_vals_type), POINTER :: section 610 INTEGER, INTENT(IN) :: start 611 612 CHARACTER(len=*), PARAMETER :: routineN = 'read_quip_section', & 613 routineP = moduleN//':'//routineN 614 615 CHARACTER(LEN=default_string_length), & 616 DIMENSION(:), POINTER :: args_str, atm_names 617 INTEGER :: is, isec, n_calc_args, n_items 618 619 CALL section_vals_get(section, n_repetition=n_items) 620 DO isec = 1, n_items 621 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 622 623 nonbonded%pot(start + isec)%pot%type = quip_type 624 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 625 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 626 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 627 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 628 CALL section_vals_val_get(section, "PARM_FILE_NAME", i_rep_section=isec, & 629 c_val=nonbonded%pot(start + isec)%pot%set(1)%quip%quip_file_name) 630 CALL section_vals_val_get(section, "INIT_ARGS", i_rep_section=isec, & 631 c_vals=args_str) 632 nonbonded%pot(start + isec)%pot%set(1)%quip%init_args = "" 633 DO is = 1, SIZE(args_str) 634 nonbonded%pot(start + isec)%pot%set(1)%quip%init_args = & 635 TRIM(nonbonded%pot(start + isec)%pot%set(1)%quip%init_args)// & 636 " "//TRIM(args_str(is)) 637 END DO ! is 638 CALL section_vals_val_get(section, "CALC_ARGS", i_rep_section=isec, & 639 n_rep_val=n_calc_args) 640 IF (n_calc_args > 0) THEN 641 CALL section_vals_val_get(section, "CALC_ARGS", i_rep_section=isec, & 642 c_vals=args_str) 643 DO is = 1, SIZE(args_str) 644 nonbonded%pot(start + isec)%pot%set(1)%quip%calc_args = & 645 TRIM(nonbonded%pot(start + isec)%pot%set(1)%quip%calc_args)// & 646 " "//TRIM(args_str(is)) 647 END DO ! is 648 END IF 649 nonbonded%pot(start + isec)%pot%rcutsq = 0.0_dp 650 END DO 651 END SUBROUTINE read_quip_section 652 653! ************************************************************************************************** 654!> \brief Reads the LJ section 655!> \param nonbonded ... 656!> \param section ... 657!> \param start ... 658!> \author teo 659! ************************************************************************************************** 660 SUBROUTINE read_lj_section(nonbonded, section, start) 661 TYPE(pair_potential_p_type), POINTER :: nonbonded 662 TYPE(section_vals_type), POINTER :: section 663 INTEGER, INTENT(IN) :: start 664 665 CHARACTER(len=*), PARAMETER :: routineN = 'read_lj_section', & 666 routineP = moduleN//':'//routineN 667 668 CHARACTER(LEN=default_string_length), & 669 DIMENSION(:), POINTER :: atm_names 670 INTEGER :: isec, n_items, n_rep 671 REAL(KIND=dp) :: epsilon, rcut, sigma 672 673 CALL section_vals_get(section, n_repetition=n_items) 674 DO isec = 1, n_items 675 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 676 CALL section_vals_val_get(section, "EPSILON", i_rep_section=isec, r_val=epsilon) 677 CALL section_vals_val_get(section, "SIGMA", i_rep_section=isec, r_val=sigma) 678 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut) 679 680 nonbonded%pot(start + isec)%pot%type = lj_charmm_type 681 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 682 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 683 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 684 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 685 nonbonded%pot(start + isec)%pot%set(1)%lj%epsilon = epsilon 686 nonbonded%pot(start + isec)%pot%set(1)%lj%sigma6 = sigma**6 687 nonbonded%pot(start + isec)%pot%set(1)%lj%sigma12 = sigma**12 688 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut 689 ! 690 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep) 691 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, & 692 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin) 693 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep) 694 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, & 695 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax) 696 END DO 697 END SUBROUTINE read_lj_section 698 699! ************************************************************************************************** 700!> \brief Reads the WILLIAMS section 701!> \param nonbonded ... 702!> \param section ... 703!> \param start ... 704!> \author teo 705! ************************************************************************************************** 706 SUBROUTINE read_wl_section(nonbonded, section, start) 707 TYPE(pair_potential_p_type), POINTER :: nonbonded 708 TYPE(section_vals_type), POINTER :: section 709 INTEGER, INTENT(IN) :: start 710 711 CHARACTER(len=*), PARAMETER :: routineN = 'read_wl_section', & 712 routineP = moduleN//':'//routineN 713 714 CHARACTER(LEN=default_string_length), & 715 DIMENSION(:), POINTER :: atm_names 716 INTEGER :: isec, n_items, n_rep 717 REAL(KIND=dp) :: a, b, c, rcut 718 719 CALL section_vals_get(section, n_repetition=n_items) 720 DO isec = 1, n_items 721 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 722 CALL section_vals_val_get(section, "A", i_rep_section=isec, r_val=a) 723 CALL section_vals_val_get(section, "B", i_rep_section=isec, r_val=b) 724 CALL section_vals_val_get(section, "C", i_rep_section=isec, r_val=c) 725 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut) 726 727 nonbonded%pot(start + isec)%pot%type = wl_type 728 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 729 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 730 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 731 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 732 nonbonded%pot(start + isec)%pot%set(1)%willis%a = a 733 nonbonded%pot(start + isec)%pot%set(1)%willis%b = b 734 nonbonded%pot(start + isec)%pot%set(1)%willis%c = c 735 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut 736 ! 737 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep) 738 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, & 739 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin) 740 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep) 741 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, & 742 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax) 743 END DO 744 END SUBROUTINE read_wl_section 745 746! ************************************************************************************************** 747!> \brief Reads the GOODWIN section 748!> \param nonbonded ... 749!> \param section ... 750!> \param start ... 751!> \author teo 752! ************************************************************************************************** 753 SUBROUTINE read_gd_section(nonbonded, section, start) 754 TYPE(pair_potential_p_type), POINTER :: nonbonded 755 TYPE(section_vals_type), POINTER :: section 756 INTEGER, INTENT(IN) :: start 757 758 CHARACTER(len=*), PARAMETER :: routineN = 'read_gd_section', & 759 routineP = moduleN//':'//routineN 760 761 CHARACTER(LEN=default_string_length), & 762 DIMENSION(:), POINTER :: atm_names 763 INTEGER :: isec, m, mc, n_items, n_rep 764 REAL(KIND=dp) :: d, dc, rcut, vr0 765 766 CALL section_vals_get(section, n_repetition=n_items) 767 DO isec = 1, n_items 768 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 769 CALL section_vals_val_get(section, "VR0", i_rep_section=isec, r_val=vr0) 770 CALL section_vals_val_get(section, "D", i_rep_section=isec, r_val=d) 771 CALL section_vals_val_get(section, "DC", i_rep_section=isec, r_val=dc) 772 CALL section_vals_val_get(section, "M", i_rep_section=isec, i_val=m) 773 CALL section_vals_val_get(section, "MC", i_rep_section=isec, i_val=mc) 774 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut) 775 776 nonbonded%pot(start + isec)%pot%type = gw_type 777 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 778 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 779 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 780 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 781 nonbonded%pot(start + isec)%pot%set(1)%goodwin%vr0 = vr0 782 nonbonded%pot(start + isec)%pot%set(1)%goodwin%d = d 783 nonbonded%pot(start + isec)%pot%set(1)%goodwin%dc = dc 784 nonbonded%pot(start + isec)%pot%set(1)%goodwin%m = m 785 nonbonded%pot(start + isec)%pot%set(1)%goodwin%mc = mc 786 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut 787 ! 788 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep) 789 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, & 790 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin) 791 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep) 792 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, & 793 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax) 794 END DO 795 END SUBROUTINE read_gd_section 796 797! ************************************************************************************************** 798!> \brief Reads the IPBV section 799!> \param nonbonded ... 800!> \param section ... 801!> \param start ... 802!> \author teo 803! ************************************************************************************************** 804 SUBROUTINE read_ipbv_section(nonbonded, section, start) 805 TYPE(pair_potential_p_type), POINTER :: nonbonded 806 TYPE(section_vals_type), POINTER :: section 807 INTEGER, INTENT(IN) :: start 808 809 CHARACTER(len=*), PARAMETER :: routineN = 'read_ipbv_section', & 810 routineP = moduleN//':'//routineN 811 812 CHARACTER(LEN=default_string_length), & 813 DIMENSION(:), POINTER :: atm_names 814 INTEGER :: isec, n_items, n_rep 815 REAL(KIND=dp) :: rcut 816 817 CALL section_vals_get(section, n_repetition=n_items) 818 DO isec = 1, n_items 819 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 820 nonbonded%pot(start + isec)%pot%type = ip_type 821 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 822 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 823 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 824 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 825 CALL set_IPBV_ff(nonbonded%pot(start + isec)%pot%at1, nonbonded%pot(start + isec)%pot%at2, & 826 nonbonded%pot(start + isec)%pot%set(1)%ipbv) 827 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut) 828 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2 829 ! 830 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep) 831 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, & 832 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin) 833 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep) 834 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, & 835 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax) 836 END DO 837 END SUBROUTINE read_ipbv_section 838 839! ************************************************************************************************** 840!> \brief Reads the BMHFT section 841!> \param nonbonded ... 842!> \param section ... 843!> \param start ... 844!> \author teo 845! ************************************************************************************************** 846 SUBROUTINE read_bmhft_section(nonbonded, section, start) 847 TYPE(pair_potential_p_type), POINTER :: nonbonded 848 TYPE(section_vals_type), POINTER :: section 849 INTEGER, INTENT(IN) :: start 850 851 CHARACTER(len=*), PARAMETER :: routineN = 'read_bmhft_section', & 852 routineP = moduleN//':'//routineN 853 854 CHARACTER(LEN=default_string_length), DIMENSION(2) :: map_atoms 855 CHARACTER(LEN=default_string_length), & 856 DIMENSION(:), POINTER :: atm_names 857 INTEGER :: i, isec, n_items, n_rep 858 REAL(KIND=dp) :: rcut 859 860 CALL section_vals_get(section, n_repetition=n_items) 861 DO isec = 1, n_items 862 CALL cite_reference(Tosi1964a) 863 CALL cite_reference(Tosi1964b) 864 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 865 nonbonded%pot(start + isec)%pot%type = ft_type 866 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 867 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 868 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 869 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 870 871 CALL section_vals_val_get(section, "A", i_rep_section=isec, n_rep_val=i) 872 IF (i == 1) THEN 873 CALL section_vals_val_get(section, "A", i_rep_section=isec, & 874 r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%a) 875 CALL section_vals_val_get(section, "B", i_rep_section=isec, & 876 r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%b) 877 CALL section_vals_val_get(section, "C", i_rep_section=isec, & 878 r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%c) 879 CALL section_vals_val_get(section, "D", i_rep_section=isec, & 880 r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%d) 881 ELSE 882 CALL section_vals_val_get(section, "MAP_ATOMS", i_rep_section=isec, c_vals=atm_names) 883 map_atoms = atm_names 884 CALL uppercase(map_atoms(1)) 885 CALL uppercase(map_atoms(2)) 886 CALL set_BMHFT_ff(map_atoms(1), map_atoms(2), nonbonded%pot(start + isec)%pot%set(1)%ft) 887 END IF 888 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut) 889 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2 890 ! 891 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep) 892 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, & 893 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin) 894 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep) 895 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, & 896 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax) 897 END DO 898 END SUBROUTINE read_bmhft_section 899 900! ************************************************************************************************** 901!> \brief Reads the BMHFTD section 902!> \param nonbonded ... 903!> \param section ... 904!> \param start ... 905!> \author Mathieu Salanne 05.2010 906! ************************************************************************************************** 907 SUBROUTINE read_bmhftd_section(nonbonded, section, start) 908 TYPE(pair_potential_p_type), POINTER :: nonbonded 909 TYPE(section_vals_type), POINTER :: section 910 INTEGER, INTENT(IN) :: start 911 912 CHARACTER(len=*), PARAMETER :: routineN = 'read_bmhftd_section', & 913 routineP = moduleN//':'//routineN 914 915 CHARACTER(LEN=default_string_length), DIMENSION(2) :: map_atoms 916 CHARACTER(LEN=default_string_length), & 917 DIMENSION(:), POINTER :: atm_names 918 INTEGER :: i, isec, n_items, n_rep 919 REAL(KIND=dp) :: rcut 920 921 CALL section_vals_get(section, n_repetition=n_items) 922 DO isec = 1, n_items 923 CALL cite_reference(Tosi1964a) 924 CALL cite_reference(Tosi1964b) 925 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 926 nonbonded%pot(start + isec)%pot%type = ftd_type 927 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 928 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 929 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 930 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 931 932 CALL section_vals_val_get(section, "A", i_rep_section=isec, n_rep_val=i) 933 IF (i == 1) THEN 934 CALL section_vals_val_get(section, "A", i_rep_section=isec, & 935 r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%a) 936 CALL section_vals_val_get(section, "B", i_rep_section=isec, & 937 r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%b) 938 CALL section_vals_val_get(section, "C", i_rep_section=isec, & 939 r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%c) 940 CALL section_vals_val_get(section, "D", i_rep_section=isec, & 941 r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%d) 942 CALL section_vals_val_get(section, "BD", i_rep_section=isec, & 943 r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%bd) 944 945 ELSE 946 CALL section_vals_val_get(section, "MAP_ATOMS", i_rep_section=isec, c_vals=atm_names) 947 map_atoms = atm_names 948 CALL uppercase(map_atoms(1)) 949 CALL uppercase(map_atoms(2)) 950 CALL set_BMHFTD_ff() 951 END IF 952 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut) 953 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2 954 ! 955 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep) 956 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, & 957 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin) 958 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep) 959 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, & 960 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax) 961 END DO 962 END SUBROUTINE read_bmhftd_section 963 964! ************************************************************************************************** 965!> \brief Reads the Buckingham 4 Ranges potential section 966!> \param nonbonded ... 967!> \param section ... 968!> \param start ... 969!> \par History 970!> MK (11.11.2010): Automatic fit of the (default) polynomial coefficients 971!> \author MI,MK 972! ************************************************************************************************** 973 SUBROUTINE read_b4_section(nonbonded, section, start) 974 975 TYPE(pair_potential_p_type), POINTER :: nonbonded 976 TYPE(section_vals_type), POINTER :: section 977 INTEGER, INTENT(IN) :: start 978 979 CHARACTER(len=*), PARAMETER :: routineN = 'read_b4_section', & 980 routineP = moduleN//':'//routineN 981 982 CHARACTER(LEN=default_string_length), & 983 DIMENSION(:), POINTER :: atm_names 984 INTEGER :: i, ir, isec, n_items, n_rep, np1, np2 985 LOGICAL :: explicit_poly1, explicit_poly2 986 REAL(KIND=dp) :: a, b, c, eval_error, r1, r2, r3, rcut 987 REAL(KIND=dp), DIMENSION(10) :: v, x 988 REAL(KIND=dp), DIMENSION(10, 10) :: p, p_inv 989 REAL(KIND=dp), DIMENSION(:), POINTER :: coeff1, coeff2, list 990 991 NULLIFY (coeff1) 992 NULLIFY (coeff2) 993 994 CALL section_vals_get(section, n_repetition=n_items) 995 996 DO isec = 1, n_items 997 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 998 CALL section_vals_val_get(section, "A", i_rep_section=isec, r_val=a) 999 CALL section_vals_val_get(section, "B", i_rep_section=isec, r_val=b) 1000 CALL section_vals_val_get(section, "C", i_rep_section=isec, r_val=c) 1001 CALL section_vals_val_get(section, "R1", i_rep_section=isec, r_val=r1) 1002 CALL section_vals_val_get(section, "R2", i_rep_section=isec, r_val=r2) 1003 CALL section_vals_val_get(section, "R3", i_rep_section=isec, r_val=r3) 1004 CALL section_vals_val_get(section, "POLY1", explicit=explicit_poly1, n_rep_val=n_rep) 1005 ! Check if polynomial coefficients were specified for range 2 and 3 explicitly 1006 IF (explicit_poly1) THEN 1007 np1 = 0 1008 DO ir = 1, n_rep 1009 NULLIFY (list) 1010 CALL section_vals_val_get(section, "POLY1", i_rep_val=ir, r_vals=list) 1011 IF (ASSOCIATED(list)) THEN 1012 CALL reallocate(coeff1, 0, np1 + SIZE(list) - 1) 1013 DO i = 1, SIZE(list) 1014 coeff1(i + np1 - 1) = list(i) 1015 END DO 1016 np1 = np1 + SIZE(list) 1017 END IF 1018 END DO 1019 END IF 1020 CALL section_vals_val_get(section, "POLY2", explicit=explicit_poly2, n_rep_val=n_rep) 1021 IF (explicit_poly2) THEN 1022 np2 = 0 1023 DO ir = 1, n_rep 1024 NULLIFY (list) 1025 CALL section_vals_val_get(section, "POLY2", i_rep_val=ir, r_vals=list) 1026 IF (ASSOCIATED(list)) THEN 1027 CALL reallocate(coeff2, 0, np2 + SIZE(list) - 1) 1028 DO i = 1, SIZE(list) 1029 coeff2(i + np2 - 1) = list(i) 1030 END DO 1031 np2 = np2 + SIZE(list) 1032 END IF 1033 END DO 1034 END IF 1035 ! Default is a 5th/3rd-order polynomial fit 1036 IF ((.NOT. explicit_poly1) .OR. (.NOT. explicit_poly2)) THEN 1037 ! Build matrix p and vector v to calculate the polynomial coefficients 1038 ! in the vector x from p*x = v 1039 p(:, :) = 0.0_dp 1040 ! Row 1: Match the 5th-order polynomial and the potential at r1 1041 p(1, 1) = 1.0_dp 1042 DO i = 2, 6 1043 p(1, i) = p(1, i - 1)*r1 1044 END DO 1045 ! Row 2: Match the first derivatives of the 5th-order polynomial and the potential at r1 1046 DO i = 2, 6 1047 p(2, i) = REAL(i - 1, KIND=dp)*p(1, i - 1) 1048 END DO 1049 ! Row 3: Match the second derivatives of the 5th-order polynomial and the potential at r1 1050 DO i = 3, 6 1051 p(3, i) = REAL(i - 1, KIND=dp)*p(2, i - 1) 1052 END DO 1053 ! Row 4: Match the 5th-order and the 3rd-order polynomials at r2 1054 p(4, 1) = 1.0_dp 1055 DO i = 2, 6 1056 p(4, i) = p(4, i - 1)*r2 1057 END DO 1058 p(4, 7) = -1.0_dp 1059 DO i = 8, 10 1060 p(4, i) = p(4, i - 1)*r2 1061 END DO 1062 ! Row 5: Match the first derivatives of the 5th-order and the 3rd-order polynomials at r2 1063 DO i = 2, 6 1064 p(5, i) = REAL(i - 1, KIND=dp)*p(4, i - 1) 1065 END DO 1066 DO i = 8, 10 1067 p(5, i) = REAL(i - 7, KIND=dp)*p(4, i - 1) 1068 END DO 1069 ! Row 6: Match the second derivatives of the 5th-order and the 3rd-order polynomials at r2 1070 DO i = 3, 6 1071 p(6, i) = REAL(i - 1, KIND=dp)*p(5, i - 1) 1072 END DO 1073 DO i = 9, 10 1074 p(6, i) = REAL(i - 7, KIND=dp)*p(5, i - 1) 1075 END DO 1076 ! Row 7: Minimum at r2, ie. the first derivative of the 3rd-order polynomial has to be zero at r2 1077 DO i = 8, 10 1078 p(7, i) = -p(5, i) 1079 END DO 1080 ! Row 8: Match the 3rd-order polynomial and the potential at r3 1081 p(8, 7) = 1.0_dp 1082 DO i = 8, 10 1083 p(8, i) = p(8, i - 1)*r3 1084 END DO 1085 ! Row 9: Match the first derivatives of the 3rd-order polynomial and the potential at r3 1086 DO i = 8, 10 1087 p(9, i) = REAL(i - 7, KIND=dp)*p(8, i - 1) 1088 END DO 1089 ! Row 10: Match the second derivatives of the 3rd-order polynomial and the potential at r3 1090 DO i = 9, 10 1091 p(10, i) = REAL(i - 7, KIND=dp)*p(9, i - 1) 1092 END DO 1093 ! Build the vector v 1094 v(1) = a*EXP(-b*r1) 1095 v(2) = -b*v(1) 1096 v(3) = -b*v(2) 1097 v(4:7) = 0.0_dp 1098 v(8) = -c/p(8, 10)**2 ! = -c/r3**6 1099 v(9) = -6.0_dp*v(8)/r3 1100 v(10) = -7.0_dp*v(9)/r3 1101 ! Calculate p_inv the inverse of the matrix p 1102 p_inv(:, :) = 0.0_dp 1103 CALL invert_matrix(p, p_inv, eval_error) 1104 IF (eval_error >= 1.0E-8_dp) & 1105 CALL cp_warn(__LOCATION__, & 1106 "The polynomial fit for the BUCK4RANGES potential is only accurate to "// & 1107 TRIM(cp_to_string(eval_error))) 1108 ! Get the 6 coefficients of the 5th-order polynomial -> x(1:6) 1109 ! and the 4 coefficients of the 3rd-order polynomial -> x(7:10) 1110 x(:) = MATMUL(p_inv(:, :), v(:)) 1111 END IF 1112 1113 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut) 1114 1115 nonbonded%pot(start + isec)%pot%type = b4_type 1116 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 1117 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 1118 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 1119 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 1120 nonbonded%pot(start + isec)%pot%set(1)%buck4r%a = a 1121 nonbonded%pot(start + isec)%pot%set(1)%buck4r%b = b 1122 nonbonded%pot(start + isec)%pot%set(1)%buck4r%c = c 1123 nonbonded%pot(start + isec)%pot%set(1)%buck4r%r1 = r1 1124 nonbonded%pot(start + isec)%pot%set(1)%buck4r%r2 = r2 1125 nonbonded%pot(start + isec)%pot%set(1)%buck4r%r3 = r3 1126 IF ((.NOT. explicit_poly1) .OR. (.NOT. explicit_poly2)) THEN 1127 nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly1 = 5 1128 nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly1(0:5) = x(1:6) 1129 nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly2 = 3 1130 nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly2(0:3) = x(7:10) 1131 ELSE 1132 nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly1 = np1 - 1 1133 CPASSERT(np1 - 1 <= 10) 1134 nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly1(0:np1 - 1) = coeff1(0:np1 - 1) 1135 nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly2 = np2 - 1 1136 CPASSERT(np2 - 1 <= 10) 1137 nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly2(0:np2 - 1) = coeff2(0:np2 - 1) 1138 END IF 1139 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut 1140 1141 IF (ASSOCIATED(coeff1)) THEN 1142 DEALLOCATE (coeff1) 1143 END IF 1144 IF (ASSOCIATED(coeff2)) THEN 1145 DEALLOCATE (coeff2) 1146 END IF 1147 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep) 1148 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, & 1149 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin) 1150 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep) 1151 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, & 1152 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax) 1153 END DO 1154 1155 END SUBROUTINE read_b4_section 1156 1157! ************************************************************************************************** 1158!> \brief Reads the GENPOT - generic potential section 1159!> \param nonbonded ... 1160!> \param section ... 1161!> \param start ... 1162!> \author Teodoro Laino - 10.2006 1163! ************************************************************************************************** 1164 SUBROUTINE read_gp_section(nonbonded, section, start) 1165 TYPE(pair_potential_p_type), POINTER :: nonbonded 1166 TYPE(section_vals_type), POINTER :: section 1167 INTEGER, INTENT(IN) :: start 1168 1169 CHARACTER(len=*), PARAMETER :: routineN = 'read_gp_section', & 1170 routineP = moduleN//':'//routineN 1171 1172 CHARACTER(LEN=default_string_length), & 1173 DIMENSION(:), POINTER :: atm_names 1174 INTEGER :: isec, n_items, n_rep 1175 REAL(KIND=dp) :: rcut 1176 1177 CALL section_vals_get(section, n_repetition=n_items) 1178 DO isec = 1, n_items 1179 NULLIFY (atm_names) 1180 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 1181 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut) 1182 nonbonded%pot(start + isec)%pot%type = gp_type 1183 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 1184 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 1185 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut 1186 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 1187 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 1188 ! Parse the genpot info 1189 CALL get_generic_info(section, "FUNCTION", nonbonded%pot(start + isec)%pot%set(1)%gp%potential, & 1190 nonbonded%pot(start + isec)%pot%set(1)%gp%parameters, & 1191 nonbonded%pot(start + isec)%pot%set(1)%gp%values, & 1192 size_variables=1, i_rep_sec=isec) 1193 nonbonded%pot(start + isec)%pot%set(1)%gp%variables = nonbonded%pot(start + isec)%pot%set(1)%gp%parameters(1) 1194 ! 1195 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep) 1196 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, & 1197 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin) 1198 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep) 1199 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, & 1200 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax) 1201 END DO 1202 END SUBROUTINE read_gp_section 1203 1204! ************************************************************************************************** 1205!> \brief Reads the tersoff section 1206!> \param nonbonded ... 1207!> \param section ... 1208!> \param start ... 1209!> \param tersoff_section ... 1210!> \author ikuo 1211! ************************************************************************************************** 1212 SUBROUTINE read_tersoff_section(nonbonded, section, start, tersoff_section) 1213 TYPE(pair_potential_p_type), POINTER :: nonbonded 1214 TYPE(section_vals_type), POINTER :: section 1215 INTEGER, INTENT(IN) :: start 1216 TYPE(section_vals_type), POINTER :: tersoff_section 1217 1218 CHARACTER(len=*), PARAMETER :: routineN = 'read_tersoff_section', & 1219 routineP = moduleN//':'//routineN 1220 1221 CHARACTER(LEN=default_string_length), & 1222 DIMENSION(:), POINTER :: atm_names 1223 INTEGER :: isec, n_items, n_rep 1224 REAL(KIND=dp) :: rcut, rcutsq 1225 1226 CALL section_vals_get(section, n_repetition=n_items) 1227 DO isec = 1, n_items 1228 CALL cite_reference(Tersoff1988) 1229 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 1230 1231 nonbonded%pot(start + isec)%pot%type = tersoff_type 1232 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 1233 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 1234 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 1235 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 1236 1237 CALL section_vals_val_get(tersoff_section, "A", i_rep_section=isec, & 1238 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%A) 1239 CALL section_vals_val_get(tersoff_section, "B", i_rep_section=isec, & 1240 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%B) 1241 CALL section_vals_val_get(tersoff_section, "lambda1", i_rep_section=isec, & 1242 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%lambda1) 1243 CALL section_vals_val_get(tersoff_section, "lambda2", i_rep_section=isec, & 1244 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%lambda2) 1245 CALL section_vals_val_get(tersoff_section, "alpha", i_rep_section=isec, & 1246 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%alpha) 1247 CALL section_vals_val_get(tersoff_section, "beta", i_rep_section=isec, & 1248 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%beta) 1249 CALL section_vals_val_get(tersoff_section, "n", i_rep_section=isec, & 1250 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%n) 1251 CALL section_vals_val_get(tersoff_section, "c", i_rep_section=isec, & 1252 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%c) 1253 CALL section_vals_val_get(tersoff_section, "d", i_rep_section=isec, & 1254 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%d) 1255 CALL section_vals_val_get(tersoff_section, "h", i_rep_section=isec, & 1256 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%h) 1257 CALL section_vals_val_get(tersoff_section, "lambda3", i_rep_section=isec, & 1258 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%lambda3) 1259 CALL section_vals_val_get(tersoff_section, "bigR", i_rep_section=isec, & 1260 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigR) 1261 CALL section_vals_val_get(tersoff_section, "bigD", i_rep_section=isec, & 1262 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigD) 1263 1264 rcutsq = (nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigR + & 1265 nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigD)**2 1266 nonbonded%pot(start + isec)%pot%set(1)%tersoff%rcutsq = rcutsq 1267 nonbonded%pot(start + isec)%pot%rcutsq = rcutsq 1268 1269 ! In case it is defined override the standard specification of RCUT 1270 CALL section_vals_val_get(tersoff_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep) 1271 IF (n_rep == 1) THEN 1272 CALL section_vals_val_get(tersoff_section, "RCUT", i_rep_section=isec, r_val=rcut) 1273 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2 1274 END IF 1275 END DO 1276 END SUBROUTINE read_tersoff_section 1277 1278! ************************************************************************************************** 1279!> \brief Reads the siepmann section 1280!> \param nonbonded ... 1281!> \param section ... 1282!> \param start ... 1283!> \param siepmann_section ... 1284!> \author Dorothea Golze 1285! ************************************************************************************************** 1286 SUBROUTINE read_siepmann_section(nonbonded, section, start, siepmann_section) 1287 TYPE(pair_potential_p_type), POINTER :: nonbonded 1288 TYPE(section_vals_type), POINTER :: section 1289 INTEGER, INTENT(IN) :: start 1290 TYPE(section_vals_type), POINTER :: siepmann_section 1291 1292 CHARACTER(len=*), PARAMETER :: routineN = 'read_siepmann_section', & 1293 routineP = moduleN//':'//routineN 1294 1295 CHARACTER(LEN=default_string_length), & 1296 DIMENSION(:), POINTER :: atm_names 1297 INTEGER :: isec, n_items, n_rep 1298 REAL(KIND=dp) :: rcut 1299 1300 CALL section_vals_get(section, n_repetition=n_items) 1301 DO isec = 1, n_items 1302 CALL cite_reference(Siepmann1995) 1303 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 1304 1305 nonbonded%pot(start + isec)%pot%type = siepmann_type 1306 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 1307 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 1308 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 1309 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 1310 1311 CALL section_vals_val_get(siepmann_section, "B", i_rep_section=isec, & 1312 r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%B) 1313 CALL section_vals_val_get(siepmann_section, "D", i_rep_section=isec, & 1314 r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%D) 1315 CALL section_vals_val_get(siepmann_section, "E", i_rep_section=isec, & 1316 r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%E) 1317 CALL section_vals_val_get(siepmann_section, "F", i_rep_section=isec, & 1318 r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%F) 1319 CALL section_vals_val_get(siepmann_section, "beta", i_rep_section=isec, & 1320 r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%beta) 1321 CALL section_vals_val_get(siepmann_section, "ALLOW_OH_FORMATION", i_rep_section=isec, & 1322 l_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%allow_oh_formation) 1323 CALL section_vals_val_get(siepmann_section, "ALLOW_H3O_FORMATION", i_rep_section=isec, & 1324 l_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%allow_h3o_formation) 1325 CALL section_vals_val_get(siepmann_section, "ALLOW_O_FORMATION", i_rep_section=isec, & 1326 l_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%allow_o_formation) 1327 1328 ! ! In case it is defined override the standard specification of RCUT 1329 CALL section_vals_val_get(siepmann_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep) 1330 IF (n_rep == 1) THEN 1331 CALL section_vals_val_get(siepmann_section, "RCUT", i_rep_section=isec, r_val=rcut) 1332 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2 1333 nonbonded%pot(start + isec)%pot%set(1)%siepmann%rcutsq = rcut**2 1334 END IF 1335 END DO 1336 END SUBROUTINE read_siepmann_section 1337 1338! ************************************************************************************************** 1339!> \brief Reads the Buckingham plus Morse potential section 1340!> \param nonbonded ... 1341!> \param section ... 1342!> \param start ... 1343!> \author MI 1344! ************************************************************************************************** 1345 SUBROUTINE read_bm_section(nonbonded, section, start) 1346 TYPE(pair_potential_p_type), POINTER :: nonbonded 1347 TYPE(section_vals_type), POINTER :: section 1348 INTEGER, INTENT(IN) :: start 1349 1350 CHARACTER(len=*), PARAMETER :: routineN = 'read_bm_section', & 1351 routineP = moduleN//':'//routineN 1352 1353 CHARACTER(LEN=default_string_length), & 1354 DIMENSION(:), POINTER :: atm_names 1355 INTEGER :: isec, n_items, n_rep 1356 REAL(KIND=dp) :: a1, a2, b1, b2, beta, c, d, f0, r0, rcut 1357 1358 CALL section_vals_get(section, n_repetition=n_items) 1359 DO isec = 1, n_items 1360 CALL cite_reference(Yamada2000) 1361 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 1362 CALL section_vals_val_get(section, "F0", i_rep_section=isec, r_val=f0) 1363 CALL section_vals_val_get(section, "A1", i_rep_section=isec, r_val=a1) 1364 CALL section_vals_val_get(section, "A2", i_rep_section=isec, r_val=a2) 1365 CALL section_vals_val_get(section, "B1", i_rep_section=isec, r_val=b1) 1366 CALL section_vals_val_get(section, "B2", i_rep_section=isec, r_val=b2) 1367 CALL section_vals_val_get(section, "C", i_rep_section=isec, r_val=c) 1368 CALL section_vals_val_get(section, "D", i_rep_section=isec, r_val=d) 1369 CALL section_vals_val_get(section, "R0", i_rep_section=isec, r_val=r0) 1370 CALL section_vals_val_get(section, "Beta", i_rep_section=isec, r_val=beta) 1371 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut) 1372 1373 nonbonded%pot(start + isec)%pot%type = bm_type 1374 nonbonded%pot(start + isec)%pot%at1 = atm_names(1) 1375 nonbonded%pot(start + isec)%pot%at2 = atm_names(2) 1376 CALL uppercase(nonbonded%pot(start + isec)%pot%at1) 1377 CALL uppercase(nonbonded%pot(start + isec)%pot%at2) 1378 nonbonded%pot(start + isec)%pot%set(1)%buckmo%f0 = f0 1379 nonbonded%pot(start + isec)%pot%set(1)%buckmo%a1 = a1 1380 nonbonded%pot(start + isec)%pot%set(1)%buckmo%a2 = a2 1381 nonbonded%pot(start + isec)%pot%set(1)%buckmo%b1 = b1 1382 nonbonded%pot(start + isec)%pot%set(1)%buckmo%b2 = b2 1383 nonbonded%pot(start + isec)%pot%set(1)%buckmo%c = c 1384 nonbonded%pot(start + isec)%pot%set(1)%buckmo%d = d 1385 nonbonded%pot(start + isec)%pot%set(1)%buckmo%r0 = r0 1386 nonbonded%pot(start + isec)%pot%set(1)%buckmo%beta = beta 1387 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut 1388 ! 1389 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep) 1390 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, & 1391 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin) 1392 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep) 1393 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, & 1394 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax) 1395 END DO 1396 END SUBROUTINE read_bm_section 1397 1398! ************************************************************************************************** 1399!> \brief Reads the CHARGE section 1400!> \param charge_atm ... 1401!> \param charge ... 1402!> \param section ... 1403!> \param start ... 1404!> \author teo 1405! ************************************************************************************************** 1406 SUBROUTINE read_chrg_section(charge_atm, charge, section, start) 1407 CHARACTER(LEN=default_string_length), & 1408 DIMENSION(:), POINTER :: charge_atm 1409 REAL(KIND=dp), DIMENSION(:), POINTER :: charge 1410 TYPE(section_vals_type), POINTER :: section 1411 INTEGER, INTENT(IN) :: start 1412 1413 CHARACTER(len=*), PARAMETER :: routineN = 'read_chrg_section', & 1414 routineP = moduleN//':'//routineN 1415 1416 CHARACTER(LEN=default_string_length) :: atm_name 1417 INTEGER :: isec, n_items 1418 1419 CALL section_vals_get(section, n_repetition=n_items) 1420 DO isec = 1, n_items 1421 CALL section_vals_val_get(section, "ATOM", i_rep_section=isec, c_val=atm_name) 1422 charge_atm(start + isec) = atm_name 1423 CALL uppercase(charge_atm(start + isec)) 1424 CALL section_vals_val_get(section, "CHARGE", i_rep_section=isec, r_val=charge(start + isec)) 1425 END DO 1426 END SUBROUTINE read_chrg_section 1427 1428! ************************************************************************************************** 1429!> \brief Reads the POLARIZABILITY section 1430!> \param apol_atm ... 1431!> \param apol ... 1432!> \param damping_list ... 1433!> \param section ... 1434!> \param start ... 1435!> \author Marcel Baer 1436! ************************************************************************************************** 1437 SUBROUTINE read_apol_section(apol_atm, apol, damping_list, section, & 1438 start) 1439 CHARACTER(LEN=default_string_length), & 1440 DIMENSION(:), POINTER :: apol_atm 1441 REAL(KIND=dp), DIMENSION(:), POINTER :: apol 1442 TYPE(damping_info_type), DIMENSION(:), POINTER :: damping_list 1443 TYPE(section_vals_type), POINTER :: section 1444 INTEGER, INTENT(IN) :: start 1445 1446 CHARACTER(len=*), PARAMETER :: routineN = 'read_apol_section', & 1447 routineP = moduleN//':'//routineN 1448 1449 CHARACTER(LEN=default_string_length) :: atm_name 1450 INTEGER :: isec, isec_damp, n_damp, n_items, & 1451 start_damp, tmp_damp 1452 TYPE(section_vals_type), POINTER :: tmp_section 1453 1454 CALL section_vals_get(section, n_repetition=n_items) 1455 NULLIFY (tmp_section) 1456 n_damp = 0 1457! *** Counts number of DIPOLE%DAMPING sections **** 1458 DO isec = 1, n_items 1459 tmp_section => section_vals_get_subs_vals(section, "DAMPING", & 1460 i_rep_section=isec) 1461 CALL section_vals_get(tmp_section, n_repetition=tmp_damp) 1462 n_damp = n_damp + tmp_damp 1463 1464 END DO 1465 1466 IF (n_damp > 0) THEN 1467 ALLOCATE (damping_list(1:n_damp)) 1468 ENDIF 1469 1470! *** Reads DIPOLE sections ***** 1471 start_damp = 0 1472 DO isec = 1, n_items 1473 CALL section_vals_val_get(section, "ATOM", i_rep_section=isec, c_val=atm_name) 1474 apol_atm(start + isec) = atm_name 1475 CALL uppercase(apol_atm(start + isec)) 1476 CALL section_vals_val_get(section, "APOL", i_rep_section=isec, r_val=apol(start + isec)) 1477 1478 tmp_section => section_vals_get_subs_vals(section, "DAMPING", & 1479 i_rep_section=isec) 1480 CALL section_vals_get(tmp_section, n_repetition=tmp_damp) 1481 DO isec_damp = 1, tmp_damp 1482 damping_list(start_damp + isec_damp)%atm_name1 = apol_atm(start + isec) 1483 CALL section_vals_val_get(tmp_section, "ATOM", i_rep_section=isec_damp, & 1484 c_val=atm_name) 1485 damping_list(start_damp + isec_damp)%atm_name2 = atm_name 1486 CALL uppercase(damping_list(start_damp + isec_damp)%atm_name2) 1487 CALL section_vals_val_get(tmp_section, "TYPE", i_rep_section=isec_damp, & 1488 c_val=atm_name) 1489 damping_list(start_damp + isec_damp)%dtype = atm_name 1490 CALL uppercase(damping_list(start_damp + isec_damp)%dtype) 1491 1492 CALL section_vals_val_get(tmp_section, "ORDER", i_rep_section=isec_damp, & 1493 i_val=damping_list(start_damp + isec_damp)%order) 1494 CALL section_vals_val_get(tmp_section, "BIJ", i_rep_section=isec_damp, & 1495 r_val=damping_list(start_damp + isec_damp)%bij) 1496 CALL section_vals_val_get(tmp_section, "CIJ", i_rep_section=isec_damp, & 1497 r_val=damping_list(start_damp + isec_damp)%cij) 1498 END DO 1499 start_damp = start_damp + tmp_damp 1500 1501 END DO 1502 1503 END SUBROUTINE read_apol_section 1504 1505! ************************************************************************************************** 1506!> \brief Reads the QUADRUPOLE POLARIZABILITY section 1507!> \param cpol_atm ... 1508!> \param cpol ... 1509!> \param section ... 1510!> \param start ... 1511!> \author Marcel Baer 1512! ************************************************************************************************** 1513 SUBROUTINE read_cpol_section(cpol_atm, cpol, section, start) 1514 CHARACTER(LEN=default_string_length), & 1515 DIMENSION(:), POINTER :: cpol_atm 1516 REAL(KIND=dp), DIMENSION(:), POINTER :: cpol 1517 TYPE(section_vals_type), POINTER :: section 1518 INTEGER, INTENT(IN) :: start 1519 1520 CHARACTER(len=*), PARAMETER :: routineN = 'read_cpol_section', & 1521 routineP = moduleN//':'//routineN 1522 1523 CHARACTER(LEN=default_string_length) :: atm_name 1524 INTEGER :: isec, n_items 1525 1526 CALL section_vals_get(section, n_repetition=n_items) 1527 DO isec = 1, n_items 1528 CALL section_vals_val_get(section, "ATOM", i_rep_section=isec, c_val=atm_name) 1529 cpol_atm(start + isec) = atm_name 1530 CALL uppercase(cpol_atm(start + isec)) 1531 CALL section_vals_val_get(section, "CPOL", i_rep_section=isec, r_val=cpol(start + isec)) 1532 END DO 1533 END SUBROUTINE read_cpol_section 1534 1535! ************************************************************************************************** 1536!> \brief Reads the SHELL section 1537!> \param shell_list ... 1538!> \param section ... 1539!> \param start ... 1540!> \author Marcella Iannuzzi 1541! ************************************************************************************************** 1542 SUBROUTINE read_shell_section(shell_list, section, start) 1543 1544 TYPE(shell_p_type), DIMENSION(:), POINTER :: shell_list 1545 TYPE(section_vals_type), POINTER :: section 1546 INTEGER, INTENT(IN) :: start 1547 1548 CHARACTER(len=*), PARAMETER :: routineN = 'read_shell_section', & 1549 routineP = moduleN//':'//routineN 1550 1551 CHARACTER(LEN=default_string_length) :: atm_name 1552 INTEGER :: i_rep, n_rep 1553 REAL(dp) :: ccharge, cutoff, k, maxdist, mfrac, & 1554 scharge 1555 1556 CALL section_vals_get(section, n_repetition=n_rep) 1557 1558 DO i_rep = 1, n_rep 1559 CALL section_vals_val_get(section, "_SECTION_PARAMETERS_", & 1560 c_val=atm_name, i_rep_section=i_rep) 1561 CALL uppercase(atm_name) 1562 shell_list(start + i_rep)%atm_name = atm_name 1563 CALL section_vals_val_get(section, "CORE_CHARGE", i_rep_section=i_rep, r_val=ccharge) 1564 shell_list(start + i_rep)%shell%charge_core = ccharge 1565 CALL section_vals_val_get(section, "SHELL_CHARGE", i_rep_section=i_rep, r_val=scharge) 1566 shell_list(start + i_rep)%shell%charge_shell = scharge 1567 CALL section_vals_val_get(section, "MASS_FRACTION", i_rep_section=i_rep, r_val=mfrac) 1568 shell_list(start + i_rep)%shell%massfrac = mfrac 1569 CALL section_vals_val_get(section, "K2_SPRING", i_rep_section=i_rep, r_val=k) 1570 IF (k < 0.0_dp) THEN 1571 CALL cp_abort(__LOCATION__, & 1572 "An invalid value was specified for the force constant k2 of the core-shell "// & 1573 "spring potential") 1574 END IF 1575 shell_list(start + i_rep)%shell%k2_spring = k 1576 CALL section_vals_val_get(section, "K4_SPRING", i_rep_section=i_rep, r_val=k) 1577 IF (k < 0.0_dp) THEN 1578 CALL cp_abort(__LOCATION__, & 1579 "An invalid value was specified for the force constant k4 of the core-shell "// & 1580 "spring potential") 1581 END IF 1582 shell_list(start + i_rep)%shell%k4_spring = k 1583 CALL section_vals_val_get(section, "MAX_DISTANCE", i_rep_section=i_rep, r_val=maxdist) 1584 shell_list(start + i_rep)%shell%max_dist = maxdist 1585 CALL section_vals_val_get(section, "SHELL_CUTOFF", i_rep_section=i_rep, r_val=cutoff) 1586 shell_list(start + i_rep)%shell%shell_cutoff = cutoff 1587 END DO 1588 1589 END SUBROUTINE read_shell_section 1590 1591! ************************************************************************************************** 1592!> \brief Reads the BONDS section 1593!> \param bond_kind ... 1594!> \param bond_a ... 1595!> \param bond_b ... 1596!> \param bond_k ... 1597!> \param bond_r0 ... 1598!> \param bond_cs ... 1599!> \param section ... 1600!> \param start ... 1601!> \author teo 1602! ************************************************************************************************** 1603 SUBROUTINE read_bonds_section(bond_kind, bond_a, bond_b, bond_k, bond_r0, bond_cs, section, start) 1604 INTEGER, DIMENSION(:), POINTER :: bond_kind 1605 CHARACTER(LEN=default_string_length), & 1606 DIMENSION(:), POINTER :: bond_a, bond_b 1607 REAL(KIND=dp), DIMENSION(:, :), POINTER :: bond_k 1608 REAL(KIND=dp), DIMENSION(:), POINTER :: bond_r0, bond_cs 1609 TYPE(section_vals_type), POINTER :: section 1610 INTEGER, INTENT(IN) :: start 1611 1612 CHARACTER(len=*), PARAMETER :: routineN = 'read_bonds_section', & 1613 routineP = moduleN//':'//routineN 1614 1615 CHARACTER(LEN=default_string_length), & 1616 DIMENSION(:), POINTER :: atm_names 1617 INTEGER :: isec, k, n_items 1618 REAL(KIND=dp), DIMENSION(:), POINTER :: Kvals 1619 1620 NULLIFY (Kvals, atm_names) 1621 CALL section_vals_get(section, n_repetition=n_items) 1622 DO isec = 1, n_items 1623 CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=bond_kind(start + isec)) 1624 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 1625 bond_a(start + isec) = atm_names(1) 1626 bond_b(start + isec) = atm_names(2) 1627 CALL uppercase(bond_a(start + isec)) 1628 CALL uppercase(bond_b(start + isec)) 1629 CALL section_vals_val_get(section, "K", i_rep_section=isec, r_vals=Kvals) 1630 CPASSERT(SIZE(Kvals) <= 3) 1631 bond_k(:, start + isec) = 0.0_dp 1632 DO k = 1, SIZE(Kvals) 1633 bond_k(k, start + isec) = Kvals(k) 1634 END DO 1635 CALL section_vals_val_get(section, "R0", i_rep_section=isec, r_val=bond_r0(start + isec)) 1636 CALL section_vals_val_get(section, "CS", i_rep_section=isec, r_val=bond_cs(start + isec)) 1637 END DO 1638 END SUBROUTINE read_bonds_section 1639 1640! ************************************************************************************************** 1641!> \brief Reads the BENDS section 1642!> \param bend_kind ... 1643!> \param bend_a ... 1644!> \param bend_b ... 1645!> \param bend_c ... 1646!> \param bend_k ... 1647!> \param bend_theta0 ... 1648!> \param bend_cb ... 1649!> \param bend_r012 ... 1650!> \param bend_r032 ... 1651!> \param bend_kbs12 ... 1652!> \param bend_kbs32 ... 1653!> \param bend_kss ... 1654!> \param bend_legendre ... 1655!> \param section ... 1656!> \param start ... 1657!> \author teo 1658! ************************************************************************************************** 1659 SUBROUTINE read_bends_section(bend_kind, bend_a, bend_b, bend_c, bend_k, bend_theta0, bend_cb, & 1660 bend_r012, bend_r032, bend_kbs12, bend_kbs32, bend_kss, bend_legendre, & 1661 section, start) 1662 INTEGER, DIMENSION(:), POINTER :: bend_kind 1663 CHARACTER(LEN=default_string_length), & 1664 DIMENSION(:), POINTER :: bend_a, bend_b, bend_c 1665 REAL(KIND=dp), DIMENSION(:), POINTER :: bend_k, bend_theta0, bend_cb, bend_r012, & 1666 bend_r032, bend_kbs12, bend_kbs32, & 1667 bend_kss 1668 TYPE(legendre_data_type), DIMENSION(:), POINTER :: bend_legendre 1669 TYPE(section_vals_type), POINTER :: section 1670 INTEGER, INTENT(IN) :: start 1671 1672 CHARACTER(len=*), PARAMETER :: routineN = 'read_bends_section', & 1673 routineP = moduleN//':'//routineN 1674 1675 CHARACTER(LEN=default_string_length), & 1676 DIMENSION(:), POINTER :: atm_names 1677 INTEGER :: isec, k, n_items, n_rep 1678 REAL(KIND=dp), DIMENSION(:), POINTER :: Kvals, r_values 1679 1680 NULLIFY (Kvals, atm_names) 1681 CALL section_vals_get(section, n_repetition=n_items) 1682 bend_legendre%order = 0 1683 DO isec = 1, n_items 1684 CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=bend_kind(start + isec)) 1685 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 1686 bend_a(start + isec) = atm_names(1) 1687 bend_b(start + isec) = atm_names(2) 1688 bend_c(start + isec) = atm_names(3) 1689 CALL uppercase(bend_a(start + isec)) 1690 CALL uppercase(bend_b(start + isec)) 1691 CALL uppercase(bend_c(start + isec)) 1692 CALL section_vals_val_get(section, "K", i_rep_section=isec, r_vals=Kvals) 1693 CPASSERT(SIZE(Kvals) == 1) 1694 bend_k(start + isec) = Kvals(1) 1695 CALL section_vals_val_get(section, "THETA0", i_rep_section=isec, r_val=bend_theta0(start + isec)) 1696 CALL section_vals_val_get(section, "CB", i_rep_section=isec, r_val=bend_cb(start + isec)) 1697 CALL section_vals_val_get(section, "R012", i_rep_section=isec, r_val=bend_r012(start + isec)) 1698 CALL section_vals_val_get(section, "R032", i_rep_section=isec, r_val=bend_r032(start + isec)) 1699 CALL section_vals_val_get(section, "KBS12", i_rep_section=isec, r_val=bend_kbs12(start + isec)) 1700 CALL section_vals_val_get(section, "KBS32", i_rep_section=isec, r_val=bend_kbs32(start + isec)) 1701 CALL section_vals_val_get(section, "KSS", i_rep_section=isec, r_val=bend_kss(start + isec)) 1702 ! get legendre based data 1703 CALL section_vals_val_get(section, "LEGENDRE", i_rep_section=isec, n_rep_val=n_rep) 1704 DO k = 1, n_rep 1705 CALL section_vals_val_get(section, "LEGENDRE", i_rep_val=k, r_vals=r_values, i_rep_section=isec) 1706 bend_legendre(start + isec)%order = SIZE(r_values) 1707 IF (ASSOCIATED(bend_legendre(start + isec)%coeffs)) THEN 1708 DEALLOCATE (bend_legendre(start + isec)%coeffs) 1709 END IF 1710 ALLOCATE (bend_legendre(start + isec)%coeffs(bend_legendre(start + isec)%order)) 1711 bend_legendre(start + isec)%coeffs = r_values 1712 END DO 1713 END DO 1714 END SUBROUTINE read_bends_section 1715 1716! ************************************************************************************************** 1717!> \brief ... 1718!> \param ub_kind ... 1719!> \param ub_a ... 1720!> \param ub_b ... 1721!> \param ub_c ... 1722!> \param ub_k ... 1723!> \param ub_r0 ... 1724!> \param section ... 1725!> \param start ... 1726! ************************************************************************************************** 1727 SUBROUTINE read_ubs_section(ub_kind, ub_a, ub_b, ub_c, ub_k, ub_r0, section, start) 1728 INTEGER, DIMENSION(:), POINTER :: ub_kind 1729 CHARACTER(LEN=default_string_length), & 1730 DIMENSION(:), POINTER :: ub_a, ub_b, ub_c 1731 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ub_k 1732 REAL(KIND=dp), DIMENSION(:), POINTER :: ub_r0 1733 TYPE(section_vals_type), POINTER :: section 1734 INTEGER, INTENT(IN) :: start 1735 1736 CHARACTER(len=*), PARAMETER :: routineN = 'read_ubs_section', & 1737 routineP = moduleN//':'//routineN 1738 1739 CHARACTER(LEN=default_string_length), & 1740 DIMENSION(:), POINTER :: atm_names 1741 INTEGER :: isec, k, n_items 1742 LOGICAL :: explicit 1743 REAL(KIND=dp), DIMENSION(:), POINTER :: Kvals 1744 TYPE(section_vals_type), POINTER :: subsection 1745 1746 NULLIFY (atm_names) 1747 CALL section_vals_get(section, n_repetition=n_items) 1748 DO isec = 1, n_items 1749 subsection => section_vals_get_subs_vals(section, "UB", i_rep_section=isec) 1750 CALL section_vals_get(subsection, explicit=explicit) 1751 IF (explicit) THEN 1752 CALL section_vals_val_get(subsection, "KIND", i_val=ub_kind(start + isec)) 1753 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 1754 ub_a(start + isec) = atm_names(1) 1755 ub_b(start + isec) = atm_names(2) 1756 ub_c(start + isec) = atm_names(3) 1757 CALL uppercase(ub_a(start + isec)) 1758 CALL uppercase(ub_b(start + isec)) 1759 CALL uppercase(ub_c(start + isec)) 1760 CALL section_vals_val_get(subsection, "K", r_vals=Kvals) 1761 CPASSERT(SIZE(Kvals) <= 3) 1762 ub_k(:, start + isec) = 0.0_dp 1763 DO k = 1, SIZE(Kvals) 1764 ub_k(k, start + isec) = Kvals(k) 1765 END DO 1766 CALL section_vals_val_get(subsection, "R0", r_val=ub_r0(start + isec)) 1767 END IF 1768 END DO 1769 END SUBROUTINE read_ubs_section 1770 1771! ************************************************************************************************** 1772!> \brief Reads the TORSIONS section 1773!> \param torsion_kind ... 1774!> \param torsion_a ... 1775!> \param torsion_b ... 1776!> \param torsion_c ... 1777!> \param torsion_d ... 1778!> \param torsion_k ... 1779!> \param torsion_phi0 ... 1780!> \param torsion_m ... 1781!> \param section ... 1782!> \param start ... 1783!> \author teo 1784! ************************************************************************************************** 1785 SUBROUTINE read_torsions_section(torsion_kind, torsion_a, torsion_b, torsion_c, torsion_d, torsion_k, & 1786 torsion_phi0, torsion_m, section, start) 1787 INTEGER, DIMENSION(:), POINTER :: torsion_kind 1788 CHARACTER(LEN=default_string_length), & 1789 DIMENSION(:), POINTER :: torsion_a, torsion_b, torsion_c, & 1790 torsion_d 1791 REAL(KIND=dp), DIMENSION(:), POINTER :: torsion_k, torsion_phi0 1792 INTEGER, DIMENSION(:), POINTER :: torsion_m 1793 TYPE(section_vals_type), POINTER :: section 1794 INTEGER, INTENT(IN) :: start 1795 1796 CHARACTER(len=*), PARAMETER :: routineN = 'read_torsions_section', & 1797 routineP = moduleN//':'//routineN 1798 1799 CHARACTER(LEN=default_string_length), & 1800 DIMENSION(:), POINTER :: atm_names 1801 INTEGER :: isec, n_items 1802 1803 NULLIFY (atm_names) 1804 CALL section_vals_get(section, n_repetition=n_items) 1805 DO isec = 1, n_items 1806 CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=torsion_kind(start + isec)) 1807 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 1808 torsion_a(start + isec) = atm_names(1) 1809 torsion_b(start + isec) = atm_names(2) 1810 torsion_c(start + isec) = atm_names(3) 1811 torsion_d(start + isec) = atm_names(4) 1812 CALL uppercase(torsion_a(start + isec)) 1813 CALL uppercase(torsion_b(start + isec)) 1814 CALL uppercase(torsion_c(start + isec)) 1815 CALL uppercase(torsion_d(start + isec)) 1816 CALL section_vals_val_get(section, "K", i_rep_section=isec, r_val=torsion_k(start + isec)) 1817 CALL section_vals_val_get(section, "PHI0", i_rep_section=isec, r_val=torsion_phi0(start + isec)) 1818 CALL section_vals_val_get(section, "M", i_rep_section=isec, i_val=torsion_m(start + isec)) 1819 ! Modify parameterisation for OPLS case 1820 IF (torsion_kind(start + isec) .EQ. do_ff_opls) THEN 1821 IF (torsion_phi0(start + isec) .NE. 0.0_dp) THEN 1822 CALL cp_warn(__LOCATION__, "PHI0 parameter was non-zero "// & 1823 "for an OPLS-type TORSION. It will be ignored.") 1824 ENDIF 1825 IF (MODULO(torsion_m(start + isec), 2) .EQ. 0) THEN 1826 ! For even M, negate the cosine using a Pi phase factor 1827 torsion_phi0(start + isec) = pi 1828 ENDIF 1829 ! the K parameter appears as K/2 in the OPLS parameterisation 1830 torsion_k(start + isec) = torsion_k(start + isec)*0.5_dp 1831 END IF 1832 END DO 1833 END SUBROUTINE read_torsions_section 1834 1835! ************************************************************************************************** 1836!> \brief Reads the IMPROPER section 1837!> \param impr_kind ... 1838!> \param impr_a ... 1839!> \param impr_b ... 1840!> \param impr_c ... 1841!> \param impr_d ... 1842!> \param impr_k ... 1843!> \param impr_phi0 ... 1844!> \param section ... 1845!> \param start ... 1846!> \author louis vanduyfhuys 1847! ************************************************************************************************** 1848 SUBROUTINE read_improper_section(impr_kind, impr_a, impr_b, impr_c, impr_d, impr_k, & 1849 impr_phi0, section, start) 1850 INTEGER, DIMENSION(:), POINTER :: impr_kind 1851 CHARACTER(LEN=default_string_length), & 1852 DIMENSION(:), POINTER :: impr_a, impr_b, impr_c, impr_d 1853 REAL(KIND=dp), DIMENSION(:), POINTER :: impr_k, impr_phi0 1854 TYPE(section_vals_type), POINTER :: section 1855 INTEGER, INTENT(IN) :: start 1856 1857 CHARACTER(len=*), PARAMETER :: routineN = 'read_improper_section', & 1858 routineP = moduleN//':'//routineN 1859 1860 CHARACTER(LEN=default_string_length), & 1861 DIMENSION(:), POINTER :: atm_names 1862 INTEGER :: isec, n_items 1863 1864 NULLIFY (atm_names) 1865 CALL section_vals_get(section, n_repetition=n_items) 1866 DO isec = 1, n_items 1867 CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=impr_kind(start + isec)) 1868 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 1869 impr_a(start + isec) = atm_names(1) 1870 impr_b(start + isec) = atm_names(2) 1871 impr_c(start + isec) = atm_names(3) 1872 impr_d(start + isec) = atm_names(4) 1873 CALL uppercase(impr_a(start + isec)) 1874 CALL uppercase(impr_b(start + isec)) 1875 CALL uppercase(impr_c(start + isec)) 1876 CALL uppercase(impr_d(start + isec)) 1877 CALL section_vals_val_get(section, "K", i_rep_section=isec, r_val=impr_k(start + isec)) 1878 CALL section_vals_val_get(section, "PHI0", i_rep_section=isec, r_val=impr_phi0(start + isec)) 1879 END DO 1880 END SUBROUTINE read_improper_section 1881 1882! ************************************************************************************************** 1883!> \brief Reads the OPBEND section 1884!> \param opbend_kind ... 1885!> \param opbend_a ... 1886!> \param opbend_b ... 1887!> \param opbend_c ... 1888!> \param opbend_d ... 1889!> \param opbend_k ... 1890!> \param opbend_phi0 ... 1891!> \param section ... 1892!> \param start ... 1893!> \author louis vanduyfhuys 1894! ************************************************************************************************** 1895 SUBROUTINE read_opbend_section(opbend_kind, opbend_a, opbend_b, opbend_c, opbend_d, opbend_k, & 1896 opbend_phi0, section, start) 1897 INTEGER, DIMENSION(:), POINTER :: opbend_kind 1898 CHARACTER(LEN=default_string_length), & 1899 DIMENSION(:), POINTER :: opbend_a, opbend_b, opbend_c, opbend_d 1900 REAL(KIND=dp), DIMENSION(:), POINTER :: opbend_k, opbend_phi0 1901 TYPE(section_vals_type), POINTER :: section 1902 INTEGER, INTENT(IN) :: start 1903 1904 CHARACTER(len=*), PARAMETER :: routineN = 'read_opbend_section', & 1905 routineP = moduleN//':'//routineN 1906 1907 CHARACTER(LEN=default_string_length), & 1908 DIMENSION(:), POINTER :: atm_names 1909 INTEGER :: isec, n_items 1910 1911 NULLIFY (atm_names) 1912 CALL section_vals_get(section, n_repetition=n_items) 1913 DO isec = 1, n_items 1914 CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=opbend_kind(start + isec)) 1915 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names) 1916 opbend_a(start + isec) = atm_names(1) 1917 opbend_b(start + isec) = atm_names(2) 1918 opbend_c(start + isec) = atm_names(3) 1919 opbend_d(start + isec) = atm_names(4) 1920 CALL uppercase(opbend_a(start + isec)) 1921 CALL uppercase(opbend_b(start + isec)) 1922 CALL uppercase(opbend_c(start + isec)) 1923 CALL uppercase(opbend_d(start + isec)) 1924 CALL section_vals_val_get(section, "K", i_rep_section=isec, r_val=opbend_k(start + isec)) 1925 CALL section_vals_val_get(section, "PHI0", i_rep_section=isec, r_val=opbend_phi0(start + isec)) 1926 END DO 1927 END SUBROUTINE read_opbend_section 1928 1929! ************************************************************************************************** 1930!> \brief Reads the force_field input section 1931!> \param ff_type ... 1932!> \param para_env ... 1933!> \param mm_section ... 1934!> \par History 1935!> JGH (30.11.2001) : moved determination of setup variables to 1936!> molecule_input 1937!> \author CJM 1938! ************************************************************************************************** 1939 SUBROUTINE read_force_field_section(ff_type, para_env, mm_section) 1940 TYPE(force_field_type), INTENT(INOUT) :: ff_type 1941 TYPE(cp_para_env_type), POINTER :: para_env 1942 TYPE(section_vals_type), POINTER :: mm_section 1943 1944 TYPE(section_vals_type), POINTER :: ff_section 1945 1946 NULLIFY (ff_section) 1947 ff_section => section_vals_get_subs_vals(mm_section, "FORCEFIELD") 1948 CALL read_force_field_section1(ff_section, mm_section, ff_type, para_env) 1949 END SUBROUTINE read_force_field_section 1950 1951! ************************************************************************************************** 1952!> \brief reads EAM potential from library 1953!> \param eam ... 1954!> \param para_env ... 1955!> \param mm_section ... 1956! ************************************************************************************************** 1957 SUBROUTINE read_eam_data(eam, para_env, mm_section) 1958 TYPE(eam_pot_type), POINTER :: eam 1959 TYPE(cp_para_env_type), POINTER :: para_env 1960 TYPE(section_vals_type), POINTER :: mm_section 1961 1962 CHARACTER(len=*), PARAMETER :: routineN = 'read_eam_data', routineP = moduleN//':'//routineN 1963 1964 INTEGER :: handle, i, iw 1965 TYPE(cp_logger_type), POINTER :: logger 1966 TYPE(cp_parser_type), POINTER :: parser 1967 1968 CALL timeset(routineN, handle) 1969 NULLIFY (parser, logger) 1970 logger => cp_get_default_logger() 1971 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", & 1972 extension=".mmLog") 1973 IF (iw > 0) WRITE (iw, *) "Reading EAM data from: ", TRIM(eam%eam_file_name) 1974 CALL parser_create(parser, TRIM(eam%eam_file_name), para_env=para_env) 1975 1976 CALL parser_get_next_line(parser, 1) 1977 IF (iw > 0) WRITE (iw, *) "Title: ", parser%input_line 1978 1979 CALL parser_get_next_line(parser, 2) 1980 READ (parser%input_line, *) eam%drar, eam%drhoar, eam%acutal, eam%npoints 1981 eam%drar = cp_unit_to_cp2k(eam%drar, "angstrom") 1982 eam%acutal = cp_unit_to_cp2k(eam%acutal, "angstrom") 1983 ! Relocating arrays with the right size 1984 CALL reallocate(eam%rho, 1, eam%npoints) 1985 CALL reallocate(eam%rhop, 1, eam%npoints) 1986 CALL reallocate(eam%rval, 1, eam%npoints) 1987 CALL reallocate(eam%rhoval, 1, eam%npoints) 1988 CALL reallocate(eam%phi, 1, eam%npoints) 1989 CALL reallocate(eam%phip, 1, eam%npoints) 1990 CALL reallocate(eam%frho, 1, eam%npoints) 1991 CALL reallocate(eam%frhop, 1, eam%npoints) 1992 ! Reading density and derivative of density (with respect to r) 1993 DO i = 1, eam%npoints 1994 CALL parser_get_next_line(parser, 1) 1995 READ (parser%input_line, *) eam%rho(i), eam%rhop(i) 1996 eam%rhop(i) = cp_unit_to_cp2k(eam%rhop(i), "angstrom^-1") 1997 eam%rval(i) = REAL(i - 1, KIND=dp)*eam%drar 1998 eam%rhoval(i) = REAL(i - 1, KIND=dp)*eam%drhoar 1999 END DO 2000 ! Reading pair potential PHI and its derivative (with respect to r) 2001 DO i = 1, eam%npoints 2002 CALL parser_get_next_line(parser, 1) 2003 READ (parser%input_line, *) eam%phi(i), eam%phip(i) 2004 eam%phi(i) = cp_unit_to_cp2k(eam%phi(i), "eV") 2005 eam%phip(i) = cp_unit_to_cp2k(eam%phip(i), "eV*angstrom^-1") 2006 END DO 2007 ! Reading embedded function and its derivative (with respect to density) 2008 DO i = 1, eam%npoints 2009 CALL parser_get_next_line(parser, 1) 2010 READ (parser%input_line, *) eam%frho(i), eam%frhop(i) 2011 eam%frho(i) = cp_unit_to_cp2k(eam%frho(i), "eV") 2012 eam%frhop(i) = cp_unit_to_cp2k(eam%frhop(i), "eV") 2013 END DO 2014 2015 IF (iw > 0) WRITE (iw, *) "Finished EAM data" 2016 CALL parser_release(parser) 2017 CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO") 2018 CALL timestop(handle) 2019 2020 END SUBROUTINE read_eam_data 2021 2022END MODULE force_fields_input 2023