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 : MODULE manybody_quip
9 :
10 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
11 : USE atomic_kind_types, ONLY: atomic_kind_type
12 : USE bibliography, ONLY: QUIP_ref, &
13 : cite_reference
14 : USE cell_types, ONLY: cell_type
15 : USE message_passing, ONLY: mp_para_env_type
16 : USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get, &
17 : fist_nonbond_env_set, &
18 : fist_nonbond_env_type, &
19 : quip_data_type
20 : USE kinds, ONLY: dp
21 : USE pair_potential_types, ONLY: pair_potential_pp_type, &
22 : pair_potential_single_type, &
23 : quip_pot_type, &
24 : quip_type
25 : USE particle_types, ONLY: particle_type
26 : USE physcon, ONLY: angstrom, &
27 : evolt
28 : #ifdef __QUIP
29 : USE quip_unified_wrapper_module, ONLY: quip_unified_wrapper
30 : #endif
31 :
32 : #include "./base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 :
36 : PRIVATE
37 :
38 : PUBLIC quip_energy_store_force_virial, quip_add_force_virial
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_quip'
41 :
42 : CONTAINS
43 :
44 : ! **************************************************************************************************
45 : !> \brief ...
46 : !> \param particle_set ...
47 : !> \param cell ...
48 : !> \param atomic_kind_set ...
49 : !> \param potparm ...
50 : !> \param fist_nonbond_env ...
51 : !> \param pot_quip ...
52 : !> \param para_env ...
53 : ! **************************************************************************************************
54 2 : SUBROUTINE quip_energy_store_force_virial(particle_set, cell, atomic_kind_set, potparm, fist_nonbond_env, &
55 : pot_quip, para_env)
56 : TYPE(particle_type), POINTER :: particle_set(:)
57 : TYPE(cell_type), POINTER :: cell
58 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
59 : TYPE(pair_potential_pp_type), POINTER :: potparm
60 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
61 : REAL(kind=dp) :: pot_quip
62 : TYPE(mp_para_env_type), OPTIONAL, &
63 : POINTER :: para_env
64 :
65 : #ifdef __QUIP
66 2 : CHARACTER(len=2), ALLOCATABLE :: elem_symbol(:)
67 : INTEGER :: i, iat, iat_use, ikind, &
68 : jkind, n_atoms, n_atoms_use, &
69 : output_unit
70 : LOGICAL :: do_parallel
71 2 : LOGICAL, ALLOCATABLE :: use_atom(:)
72 : REAL(kind=dp) :: lattice(3, 3), virial(3, 3)
73 2 : REAL(kind=dp), ALLOCATABLE :: force(:, :), pos(:, :)
74 : TYPE(pair_potential_single_type), &
75 : POINTER :: pot
76 : TYPE(quip_data_type), POINTER :: quip_data
77 : TYPE(quip_pot_type), POINTER :: quip
78 :
79 : #endif
80 : #ifndef __QUIP
81 : MARK_USED(particle_set)
82 : MARK_USED(cell)
83 : MARK_USED(atomic_kind_set)
84 : MARK_USED(potparm)
85 : MARK_USED(fist_nonbond_env)
86 : MARK_USED(pot_quip)
87 : MARK_USED(para_env)
88 : CALL cp_abort(__LOCATION__, "In order to use QUIP you need to download "// &
89 : "and install the libAtoms/QUIP library (check CP2K manual for details)")
90 : #else
91 2 : n_atoms = SIZE(particle_set)
92 6 : ALLOCATE (use_atom(n_atoms))
93 8 : use_atom = .FALSE.
94 :
95 2 : NULLIFY (quip)
96 :
97 4 : DO ikind = 1, SIZE(atomic_kind_set)
98 6 : DO jkind = 1, SIZE(atomic_kind_set)
99 2 : pot => potparm%pot(ikind, jkind)%pot
100 6 : DO i = 1, SIZE(pot%type)
101 2 : IF (pot%type(i) /= quip_type) CYCLE
102 2 : IF (.NOT. ASSOCIATED(quip)) quip => pot%set(i)%quip
103 10 : DO iat = 1, n_atoms
104 6 : IF (particle_set(iat)%atomic_kind%kind_number == ikind .OR. &
105 8 : particle_set(iat)%atomic_kind%kind_number == jkind) use_atom(iat) = .TRUE.
106 : END DO ! iat
107 : END DO ! i
108 : END DO ! jkind
109 : END DO ! ikind
110 8 : n_atoms_use = COUNT(use_atom)
111 10 : ALLOCATE (pos(3, n_atoms_use), force(3, n_atoms_use), elem_symbol(n_atoms_use))
112 :
113 8 : iat_use = 0
114 8 : DO iat = 1, n_atoms
115 6 : IF (.NOT. use_atom(iat)) CYCLE
116 6 : iat_use = iat_use + 1
117 24 : pos(1:3, iat_use) = particle_set(iat)%r*angstrom
118 8 : elem_symbol(iat_use) = particle_set(iat)%atomic_kind%element_symbol
119 : END DO
120 2 : IF (iat_use > 0) CALL cite_reference(QUIP_ref)
121 2 : output_unit = cp_logger_get_default_io_unit()
122 26 : lattice = cell%hmat*angstrom
123 2 : do_parallel = .FALSE.
124 2 : IF (PRESENT(para_env)) THEN
125 2 : do_parallel = para_env%num_pe > 1
126 : END IF
127 2 : IF (do_parallel) THEN
128 : CALL quip_unified_wrapper( &
129 : N=n_atoms_use, pos=pos, lattice=lattice, symbol=elem_symbol, &
130 : quip_param_file=TRIM(quip%quip_file_name), &
131 : quip_param_file_len=LEN_TRIM(quip%quip_file_name), &
132 : init_args_str=TRIM(quip%init_args), &
133 : init_args_str_len=LEN_TRIM(quip%init_args), &
134 : calc_args_str=TRIM(quip%calc_args), &
135 : calc_args_str_len=LEN_TRIM(quip%calc_args), &
136 : energy=pot_quip, force=force, virial=virial, &
137 2 : output_unit=output_unit, mpi_communicator=para_env%get_handle())
138 : ELSE
139 : CALL quip_unified_wrapper( &
140 : N=n_atoms_use, pos=pos, lattice=lattice, symbol=elem_symbol, &
141 : quip_param_file=TRIM(quip%quip_file_name), &
142 : quip_param_file_len=LEN_TRIM(quip%quip_file_name), &
143 : init_args_str=TRIM(quip%init_args), &
144 : init_args_str_len=LEN_TRIM(quip%init_args), &
145 : calc_args_str=TRIM(quip%calc_args), &
146 : calc_args_str_len=LEN_TRIM(quip%calc_args), &
147 0 : energy=pot_quip, force=force, virial=virial, output_unit=output_unit)
148 : END IF
149 : ! convert units
150 2 : pot_quip = pot_quip/evolt
151 26 : force = force/(evolt/angstrom)
152 26 : virial = virial/evolt
153 : ! account for double counting from multiple MPI processes
154 2 : IF (PRESENT(para_env)) pot_quip = pot_quip/REAL(para_env%num_pe, dp)
155 26 : IF (PRESENT(para_env)) force = force/REAL(para_env%num_pe, dp)
156 26 : IF (PRESENT(para_env)) virial = virial/REAL(para_env%num_pe, dp)
157 : ! get quip_data to save force, virial info
158 2 : CALL fist_nonbond_env_get(fist_nonbond_env, quip_data=quip_data)
159 2 : IF (.NOT. ASSOCIATED(quip_data)) THEN
160 26 : ALLOCATE (quip_data)
161 2 : CALL fist_nonbond_env_set(fist_nonbond_env, quip_data=quip_data)
162 2 : NULLIFY (quip_data%use_indices, quip_data%force)
163 : END IF
164 2 : IF (ASSOCIATED(quip_data%force)) THEN
165 0 : IF (SIZE(quip_data%force, 2) /= n_atoms_use) THEN
166 0 : DEALLOCATE (quip_data%force, quip_data%use_indices)
167 : END IF
168 : END IF
169 2 : IF (.NOT. ASSOCIATED(quip_data%force)) THEN
170 4 : ALLOCATE (quip_data%force(3, n_atoms_use))
171 6 : ALLOCATE (quip_data%use_indices(n_atoms_use))
172 : END IF
173 : ! save force, virial info
174 : iat_use = 0
175 8 : DO iat = 1, n_atoms
176 8 : IF (use_atom(iat)) THEN
177 6 : iat_use = iat_use + 1
178 6 : quip_data%use_indices(iat_use) = iat
179 : END IF
180 : END DO
181 26 : quip_data%force = force
182 26 : quip_data%virial = virial
183 :
184 2 : DEALLOCATE (use_atom, pos, force, elem_symbol)
185 : #endif
186 2 : END SUBROUTINE quip_energy_store_force_virial
187 :
188 : ! **************************************************************************************************
189 : !> \brief ...
190 : !> \param fist_nonbond_env ...
191 : !> \param force ...
192 : !> \param virial ...
193 : ! **************************************************************************************************
194 9330 : SUBROUTINE quip_add_force_virial(fist_nonbond_env, force, virial)
195 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
196 : REAL(KIND=dp) :: force(:, :), virial(3, 3)
197 :
198 : #ifdef __QUIP
199 : INTEGER :: iat, iat_use
200 : TYPE(quip_data_type), POINTER :: quip_data
201 : #endif
202 :
203 : #ifndef __QUIP
204 : MARK_USED(fist_nonbond_env)
205 : MARK_USED(force)
206 : MARK_USED(virial)
207 : RETURN
208 : #else
209 9330 : CALL fist_nonbond_env_get(fist_nonbond_env, quip_data=quip_data)
210 9330 : IF (.NOT. ASSOCIATED(quip_data)) RETURN
211 :
212 0 : DO iat_use = 1, SIZE(quip_data%use_indices)
213 0 : iat = quip_data%use_indices(iat_use)
214 0 : CPASSERT(iat >= 1 .AND. iat <= SIZE(force, 2))
215 0 : force(1:3, iat) = force(1:3, iat) + quip_data%force(1:3, iat_use)
216 : END DO
217 0 : virial = virial + quip_data%virial
218 : #endif
219 : END SUBROUTINE quip_add_force_virial
220 :
221 : END MODULE manybody_quip
|
