Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 fparser, ONLY: evalf,&
23 : evalfd,&
24 : finalizef,&
25 : initf,&
26 : parsef
27 : USE generic_function_input, ONLY: parse_function_section
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 102275 : 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 102275 : ALLOCATABLE, DIMENSION(:) :: my_par
63 : INTEGER :: a_var, handle, i, iatom, j, k, n_var, &
64 : natom, rep
65 102275 : INTEGER, DIMENSION(:), POINTER :: iatms, nparticle
66 : LOGICAL :: useall
67 : REAL(KIND=dp) :: dedf, dx, energy, err, lerr
68 102275 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: 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 102275 : useall = .FALSE.
75 102275 : CALL timeset(routineN, handle)
76 102275 : NULLIFY (logger, subsys, particles, ext_pot_section, nparticle)
77 : ext_pot_section => section_vals_get_subs_vals(force_env%force_env_section, &
78 102275 : "EXTERNAL_POTENTIAL")
79 102275 : CALL section_vals_get(ext_pot_section, n_repetition=n_var)
80 102669 : DO rep = 1, n_var
81 394 : natom = 0
82 394 : logger => cp_get_default_logger()
83 : CALL parse_function_section(ext_pot_section, "FUNCTION", coupling_function, &
84 : my_par, my_val, i_rep_sec=rep, &
85 1576 : input_variables=["X", "Y", "Z"], dx=dx, error_limit=lerr)
86 394 : CALL initf(1)
87 394 : CALL parsef(1, TRIM(coupling_function), my_par)
88 :
89 : ! Apply potential on all atoms, computing energy and forces
90 394 : NULLIFY (particles, subsys)
91 394 : CALL force_env_get(force_env, subsys=subsys)
92 394 : CALL cp_subsys_get(subsys, particles=particles)
93 394 : CALL force_env_get(force_env, additional_potential=energy)
94 394 : CALL section_vals_val_get(ext_pot_section, "ATOMS_LIST", n_rep_val=a_var, i_rep_section=rep)
95 504 : DO k = 1, a_var
96 110 : CALL section_vals_val_get(ext_pot_section, "ATOMS_LIST", i_rep_val=k, i_vals=iatms, i_rep_section=rep)
97 110 : CALL reallocate(nparticle, 1, natom + SIZE(iatms))
98 484 : nparticle(natom + 1:natom + SIZE(iatms)) = iatms
99 504 : natom = natom + SIZE(iatms)
100 : END DO
101 394 : IF (a_var == 0) THEN
102 284 : natom = particles%n_els
103 284 : useall = .TRUE.
104 : END IF
105 854 : DO i = 1, natom
106 460 : IF (useall) THEN
107 : iatom = i
108 : ELSE
109 132 : iatom = nparticle(i)
110 : END IF
111 460 : my_val(1) = particles%els(iatom)%r(1)
112 460 : my_val(2) = particles%els(iatom)%r(2)
113 460 : my_val(3) = particles%els(iatom)%r(3)
114 :
115 460 : energy = energy + evalf(1, my_val)
116 2234 : DO j = 1, 3
117 1380 : dedf = evalfd(1, j, my_val, dx, err)
118 1380 : IF (ABS(err) > lerr) THEN
119 110 : WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
120 110 : WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
121 110 : CALL compress(this_error, .TRUE.)
122 110 : CALL compress(def_error, .TRUE.)
123 : CALL cp_warn(__LOCATION__, &
124 : 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
125 : ' Error '//TRIM(this_error)//' in computing numerical derivatives larger than'// &
126 110 : TRIM(def_error)//' .')
127 : END IF
128 1840 : particles%els(iatom)%f(j) = particles%els(iatom)%f(j) - dedf
129 : END DO
130 : END DO
131 394 : CALL force_env_set(force_env, additional_potential=energy)
132 394 : DEALLOCATE (my_par)
133 394 : DEALLOCATE (my_val)
134 394 : IF (a_var /= 0) THEN
135 110 : DEALLOCATE (nparticle)
136 : END IF
137 103063 : CALL finalizef()
138 : END DO
139 102275 : CALL timestop(handle)
140 204550 : END SUBROUTINE add_external_potential
141 :
142 : END MODULE external_potential_methods
|
