Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Methods to include the effect of an external potential during an MD
10 : !> or energy calculation
11 : !> \author Teodoro Laino (03.2008) [tlaino]
12 : ! **************************************************************************************************
13 : MODULE external_potential_methods
14 : USE cp_log_handling, ONLY: cp_get_default_logger,&
15 : cp_logger_type,&
16 : cp_to_string
17 : USE cp_subsys_types, ONLY: cp_subsys_get,&
18 : cp_subsys_type
19 : USE force_env_types, ONLY: force_env_get,&
20 : force_env_set,&
21 : force_env_type
22 : USE force_fields_util, ONLY: get_generic_info
23 : USE fparser, ONLY: evalf,&
24 : evalfd,&
25 : finalizef,&
26 : initf,&
27 : parsef
28 : USE input_section_types, ONLY: section_vals_get,&
29 : section_vals_get_subs_vals,&
30 : section_vals_type,&
31 : section_vals_val_get
32 : USE kinds, ONLY: default_path_length,&
33 : default_string_length,&
34 : dp
35 : USE memory_utilities, ONLY: reallocate
36 : USE particle_list_types, ONLY: particle_list_type
37 : USE string_utilities, ONLY: compress
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'external_potential_methods'
44 : PUBLIC :: add_external_potential
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief ...
50 : !> \param force_env ...
51 : !> \date 03.2008
52 : !> \author Teodoro Laino - University of Zurich [tlaino]
53 : ! **************************************************************************************************
54 196024 : SUBROUTINE add_external_potential(force_env)
55 : TYPE(force_env_type), POINTER :: force_env
56 :
57 : CHARACTER(len=*), PARAMETER :: routineN = 'add_external_potential'
58 :
59 : CHARACTER(LEN=default_path_length) :: coupling_function
60 : CHARACTER(LEN=default_string_length) :: def_error, this_error
61 : CHARACTER(LEN=default_string_length), &
62 98012 : DIMENSION(:), POINTER :: my_par
63 : INTEGER :: a_var, handle, i, iatom, j, k, n_var, &
64 : natom, rep
65 98012 : INTEGER, DIMENSION(:), POINTER :: iatms, nparticle
66 : LOGICAL :: useall
67 : REAL(KIND=dp) :: dedf, dx, energy, err, lerr
68 98012 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_val
69 : TYPE(cp_logger_type), POINTER :: logger
70 : TYPE(cp_subsys_type), POINTER :: subsys
71 : TYPE(particle_list_type), POINTER :: particles
72 : TYPE(section_vals_type), POINTER :: ext_pot_section
73 :
74 98012 : useall = .FALSE.
75 98012 : CALL timeset(routineN, handle)
76 98012 : NULLIFY (my_par, my_val, logger, subsys, particles, ext_pot_section, nparticle)
77 : ext_pot_section => section_vals_get_subs_vals(force_env%force_env_section, &
78 98012 : "EXTERNAL_POTENTIAL")
79 98012 : CALL section_vals_get(ext_pot_section, n_repetition=n_var)
80 98406 : DO rep = 1, n_var
81 394 : natom = 0
82 394 : logger => cp_get_default_logger()
83 394 : CALL section_vals_val_get(ext_pot_section, "DX", r_val=dx, i_rep_section=rep)
84 394 : CALL section_vals_val_get(ext_pot_section, "ERROR_LIMIT", r_val=lerr, i_rep_section=rep)
85 : CALL get_generic_info(ext_pot_section, "FUNCTION", coupling_function, my_par, my_val, &
86 1576 : input_variables=(/"X", "Y", "Z"/), i_rep_sec=rep)
87 394 : CALL initf(1)
88 394 : CALL parsef(1, TRIM(coupling_function), my_par)
89 :
90 : ! Apply potential on all atoms, computing energy and forces
91 394 : NULLIFY (particles, subsys)
92 394 : CALL force_env_get(force_env, subsys=subsys)
93 394 : CALL cp_subsys_get(subsys, particles=particles)
94 394 : CALL force_env_get(force_env, additional_potential=energy)
95 394 : CALL section_vals_val_get(ext_pot_section, "ATOMS_LIST", n_rep_val=a_var, i_rep_section=rep)
96 504 : DO k = 1, a_var
97 110 : CALL section_vals_val_get(ext_pot_section, "ATOMS_LIST", i_rep_val=k, i_vals=iatms, i_rep_section=rep)
98 110 : CALL reallocate(nparticle, 1, natom + SIZE(iatms))
99 484 : nparticle(natom + 1:natom + SIZE(iatms)) = iatms
100 504 : natom = natom + SIZE(iatms)
101 : END DO
102 394 : IF (a_var == 0) THEN
103 284 : natom = particles%n_els
104 284 : useall = .TRUE.
105 : END IF
106 854 : DO i = 1, natom
107 460 : IF (useall) THEN
108 : iatom = i
109 : ELSE
110 132 : iatom = nparticle(i)
111 : END IF
112 460 : my_val(1) = particles%els(iatom)%r(1)
113 460 : my_val(2) = particles%els(iatom)%r(2)
114 460 : my_val(3) = particles%els(iatom)%r(3)
115 :
116 460 : energy = energy + evalf(1, my_val)
117 2234 : DO j = 1, 3
118 1380 : dedf = evalfd(1, j, my_val, dx, err)
119 1380 : IF (ABS(err) > lerr) THEN
120 110 : WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
121 110 : WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
122 110 : CALL compress(this_error, .TRUE.)
123 110 : CALL compress(def_error, .TRUE.)
124 : CALL cp_warn(__LOCATION__, &
125 : 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
126 : ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'// &
127 110 : TRIM(def_error)//' .')
128 : END IF
129 1840 : particles%els(iatom)%f(j) = particles%els(iatom)%f(j) - dedf
130 : END DO
131 : END DO
132 394 : CALL force_env_set(force_env, additional_potential=energy)
133 394 : DEALLOCATE (my_par)
134 394 : DEALLOCATE (my_val)
135 394 : IF (a_var /= 0) THEN
136 110 : DEALLOCATE (nparticle)
137 : END IF
138 99194 : CALL finalizef()
139 : END DO
140 98012 : CALL timestop(handle)
141 98012 : END SUBROUTINE add_external_potential
142 :
143 : END MODULE external_potential_methods
|
