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 Define the data structure for the molecule information.
10 : !> \par History
11 : !> JGH (22.05.2004) add last_atom information
12 : !> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
13 : !> (patch by Marcel Baer)
14 : !> \author Matthias Krack (29.08.2003)
15 : ! **************************************************************************************************
16 : MODULE molecule_types
17 :
18 : USE colvar_types, ONLY: colvar_counters,&
19 : colvar_release,&
20 : colvar_type
21 : USE kinds, ONLY: dp
22 : USE molecule_kind_types, ONLY: colvar_constraint_type,&
23 : fixd_constraint_type,&
24 : g3x3_constraint_type,&
25 : g4x6_constraint_type,&
26 : get_molecule_kind,&
27 : molecule_kind_type,&
28 : vsite_constraint_type
29 : #include "../base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 :
33 : PRIVATE
34 :
35 : ! Global parameters (in this module)
36 :
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecule_types'
38 :
39 : ! Molecular constraint types
40 : TYPE local_colvar_constraint_type
41 : TYPE(colvar_type), POINTER :: colvar => NULL(), &
42 : colvar_old => NULL()
43 : REAL(KIND=dp) :: lambda = 0.0_dp, &
44 : sigma = 0.0_dp
45 : LOGICAL :: init = .FALSE.
46 : END TYPE local_colvar_constraint_type
47 :
48 : TYPE local_g3x3_constraint_type
49 : LOGICAL :: init = .FALSE.
50 : REAL(KIND=dp) :: scale = 0.0_dp, &
51 : imass1 = 0.0_dp, &
52 : imass2 = 0.0_dp, &
53 : imass3 = 0.0_dp, &
54 : scale_old = 0.0_dp
55 : REAL(KIND=dp), DIMENSION(3) :: fa = 0.0_dp, &
56 : fb = 0.0_dp, &
57 : fc = 0.0_dp, &
58 : f_roll1 = 0.0_dp, &
59 : f_roll2 = 0.0_dp, &
60 : f_roll3 = 0.0_dp, &
61 : ra_old = 0.0_dp, &
62 : rb_old = 0.0_dp, &
63 : rc_old = 0.0_dp, &
64 : r0_12 = 0.0_dp, &
65 : r0_13 = 0.0_dp, &
66 : r0_23 = 0.0_dp, &
67 : va = 0.0_dp, &
68 : vb = 0.0_dp, &
69 : vc = 0.0_dp, &
70 : del_lambda = 0.0_dp, &
71 : lambda = 0.0_dp, &
72 : lambda_old = 0.0_dp
73 : REAL(KIND=dp), DIMENSION(3, 3) :: amat = 0.0_dp
74 : END TYPE local_g3x3_constraint_type
75 :
76 : TYPE local_g4x6_constraint_type
77 : LOGICAL :: init = .FALSE.
78 : REAL(KIND=dp) :: scale = 0.0_dp, &
79 : scale_old = 0.0_dp, &
80 : imass1 = 0.0_dp, &
81 : imass2 = 0.0_dp, &
82 : imass3 = 0.0_dp, &
83 : imass4 = 0.0_dp
84 : REAL(KIND=dp), DIMENSION(3) :: fa = 0.0_dp, &
85 : fb = 0.0_dp, &
86 : fc = 0.0_dp, &
87 : fd = 0.0_dp, &
88 : fe = 0.0_dp, &
89 : ff = 0.0_dp, &
90 : f_roll1 = 0.0_dp, &
91 : f_roll2 = 0.0_dp, &
92 : f_roll3 = 0.0_dp, &
93 : f_roll4 = 0.0_dp, &
94 : f_roll5 = 0.0_dp, &
95 : f_roll6 = 0.0_dp, &
96 : ra_old = 0.0_dp, &
97 : rb_old = 0.0_dp, &
98 : rc_old = 0.0_dp, &
99 : rd_old = 0.0_dp, &
100 : re_old = 0.0_dp, &
101 : rf_old = 0.0_dp, &
102 : va = 0.0_dp, &
103 : vb = 0.0_dp, &
104 : vc = 0.0_dp, &
105 : vd = 0.0_dp, &
106 : ve = 0.0_dp, &
107 : vf = 0.0_dp, &
108 : r0_12 = 0.0_dp, &
109 : r0_13 = 0.0_dp, &
110 : r0_14 = 0.0_dp, &
111 : r0_23 = 0.0_dp, &
112 : r0_24 = 0.0_dp, &
113 : r0_34 = 0.0_dp
114 : REAL(KIND=dp), DIMENSION(6) :: del_lambda = 0.0_dp, &
115 : lambda = 0.0_dp, &
116 : lambda_old = 0.0_dp
117 : REAL(KIND=dp), DIMENSION(6, 6) :: amat = 0.0_dp
118 : END TYPE local_g4x6_constraint_type
119 :
120 : TYPE local_states_type
121 : INTEGER :: nstates = 0 ! Kohn-Sham states for molecule
122 : INTEGER, DIMENSION(:), POINTER :: states => NULL() ! indices of Kohn-Sham states for molecule
123 : END TYPE local_states_type
124 :
125 : TYPE local_constraint_type
126 : TYPE(local_colvar_constraint_type), &
127 : DIMENSION(:), POINTER :: lcolv => NULL()
128 : TYPE(local_g3x3_constraint_type), DIMENSION(:), &
129 : POINTER :: lg3x3 => NULL()
130 : TYPE(local_g4x6_constraint_type), DIMENSION(:), &
131 : POINTER :: lg4x6 => NULL()
132 : END TYPE local_constraint_type
133 :
134 : TYPE global_constraint_type
135 : TYPE(colvar_counters) :: ncolv = colvar_counters()
136 : LOGICAL :: joint_colvar_checks = .FALSE.
137 : INTEGER :: ntot = 0, &
138 : nrestraint = 0, &
139 : ng3x3 = 0, &
140 : ng3x3_restraint = 0, &
141 : ng4x6 = 0, &
142 : ng4x6_restraint = 0, &
143 : nvsite = 0, &
144 : nvsite_restraint = 0
145 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list => NULL()
146 : TYPE(colvar_constraint_type), DIMENSION(:), &
147 : POINTER :: colv_list => NULL()
148 : TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list => NULL()
149 : TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list => NULL()
150 : TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list => NULL()
151 : TYPE(local_colvar_constraint_type), &
152 : DIMENSION(:), POINTER :: lcolv => NULL()
153 : TYPE(local_g3x3_constraint_type), DIMENSION(:), &
154 : POINTER :: lg3x3 => NULL()
155 : TYPE(local_g4x6_constraint_type), DIMENSION(:), &
156 : POINTER :: lg4x6 => NULL()
157 : END TYPE global_constraint_type
158 :
159 : ! Molecule type
160 : TYPE molecule_type
161 : TYPE(molecule_kind_type), POINTER :: molecule_kind => NULL() ! pointer to molecule kind information
162 : TYPE(local_states_type), DIMENSION(:), POINTER :: lmi => NULL() ! local (spin)-states information
163 : TYPE(local_constraint_type), POINTER :: lci => NULL() ! local molecule constraint info
164 : INTEGER :: first_atom = 0 ! global index of first atom in molecule
165 : INTEGER :: last_atom = 0 ! global index of last atom in molecule
166 : INTEGER :: first_shell = 0 ! global index of first shell atom in molecule
167 : INTEGER :: last_shell = 0 ! global index of last shell atom in molecule
168 : END TYPE molecule_type
169 :
170 : ! Public data types
171 :
172 : PUBLIC :: local_colvar_constraint_type, &
173 : local_g3x3_constraint_type, &
174 : local_g4x6_constraint_type, &
175 : local_constraint_type, &
176 : local_states_type, &
177 : global_constraint_type, &
178 : molecule_type
179 :
180 : ! Public subroutines
181 :
182 : PUBLIC :: deallocate_global_constraint, &
183 : allocate_molecule_set, &
184 : deallocate_molecule_set, &
185 : get_molecule, &
186 : set_molecule, &
187 : set_molecule_set, &
188 : molecule_of_atom, &
189 : get_molecule_set_info, &
190 : get_domain_set_info
191 :
192 : CONTAINS
193 :
194 : ! **************************************************************************************************
195 : !> \brief Deallocate a global constraint.
196 : !> \param gci ...
197 : !> \par History
198 : !> 07.2003 created [fawzi]
199 : !> 01.2014 moved from cp_subsys_release() into separate routine.
200 : !> \author Ole Schuett
201 : ! **************************************************************************************************
202 12138 : SUBROUTINE deallocate_global_constraint(gci)
203 : TYPE(global_constraint_type), POINTER :: gci
204 :
205 : INTEGER :: i
206 :
207 12138 : IF (ASSOCIATED(gci)) THEN
208 : ! List of constraints
209 11586 : IF (ASSOCIATED(gci%colv_list)) THEN
210 130 : DO i = 1, SIZE(gci%colv_list)
211 130 : DEALLOCATE (gci%colv_list(i)%i_atoms)
212 : END DO
213 54 : DEALLOCATE (gci%colv_list)
214 : END IF
215 :
216 11586 : IF (ASSOCIATED(gci%g3x3_list)) THEN
217 4 : DEALLOCATE (gci%g3x3_list)
218 : END IF
219 :
220 11586 : IF (ASSOCIATED(gci%g4x6_list)) THEN
221 4 : DEALLOCATE (gci%g4x6_list)
222 : END IF
223 :
224 : ! Local information
225 11586 : IF (ASSOCIATED(gci%lcolv)) THEN
226 130 : DO i = 1, SIZE(gci%lcolv)
227 76 : CALL colvar_release(gci%lcolv(i)%colvar)
228 130 : CALL colvar_release(gci%lcolv(i)%colvar_old)
229 : END DO
230 54 : DEALLOCATE (gci%lcolv)
231 : END IF
232 :
233 11586 : IF (ASSOCIATED(gci%lg3x3)) THEN
234 4 : DEALLOCATE (gci%lg3x3)
235 : END IF
236 :
237 11586 : IF (ASSOCIATED(gci%lg4x6)) THEN
238 4 : DEALLOCATE (gci%lg4x6)
239 : END IF
240 :
241 11586 : IF (ASSOCIATED(gci%fixd_list)) THEN
242 2 : DEALLOCATE (gci%fixd_list)
243 : END IF
244 :
245 11586 : DEALLOCATE (gci)
246 : END IF
247 12138 : END SUBROUTINE deallocate_global_constraint
248 :
249 : ! **************************************************************************************************
250 : !> \brief Allocate a molecule set.
251 : !> \param molecule_set ...
252 : !> \param nmolecule ...
253 : !> \date 29.08.2003
254 : !> \author Matthias Krack
255 : !> \version 1.0
256 : ! **************************************************************************************************
257 12139 : SUBROUTINE allocate_molecule_set(molecule_set, nmolecule)
258 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
259 : INTEGER, INTENT(IN) :: nmolecule
260 :
261 12139 : IF (ASSOCIATED(molecule_set)) CALL deallocate_molecule_set(molecule_set)
262 :
263 356103 : ALLOCATE (molecule_set(nmolecule))
264 :
265 12139 : END SUBROUTINE allocate_molecule_set
266 :
267 : ! **************************************************************************************************
268 : !> \brief Deallocate a molecule set.
269 : !> \param molecule_set ...
270 : !> \date 29.08.2003
271 : !> \author Matthias Krack
272 : !> \version 1.0
273 : ! **************************************************************************************************
274 12138 : SUBROUTINE deallocate_molecule_set(molecule_set)
275 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
276 :
277 : INTEGER :: imolecule, j
278 :
279 12138 : IF (ASSOCIATED(molecule_set)) THEN
280 :
281 331822 : DO imolecule = 1, SIZE(molecule_set)
282 319684 : IF (ASSOCIATED(molecule_set(imolecule)%lmi)) THEN
283 70 : DO j = 1, SIZE(molecule_set(imolecule)%lmi)
284 70 : IF (ASSOCIATED(molecule_set(imolecule)%lmi(j)%states)) THEN
285 40 : DEALLOCATE (molecule_set(imolecule)%lmi(j)%states)
286 : END IF
287 : END DO
288 30 : DEALLOCATE (molecule_set(imolecule)%lmi)
289 : END IF
290 331822 : IF (ASSOCIATED(molecule_set(imolecule)%lci)) THEN
291 43640 : IF (ASSOCIATED(molecule_set(imolecule)%lci%lcolv)) THEN
292 4366 : DO j = 1, SIZE(molecule_set(imolecule)%lci%lcolv)
293 2228 : CALL colvar_release(molecule_set(imolecule)%lci%lcolv(j)%colvar)
294 4366 : CALL colvar_release(molecule_set(imolecule)%lci%lcolv(j)%colvar_old)
295 : END DO
296 2138 : DEALLOCATE (molecule_set(imolecule)%lci%lcolv)
297 : END IF
298 43640 : IF (ASSOCIATED(molecule_set(imolecule)%lci%lg3x3)) THEN
299 36290 : DEALLOCATE (molecule_set(imolecule)%lci%lg3x3)
300 : END IF
301 43640 : IF (ASSOCIATED(molecule_set(imolecule)%lci%lg4x6)) THEN
302 650 : DEALLOCATE (molecule_set(imolecule)%lci%lg4x6)
303 : END IF
304 43640 : DEALLOCATE (molecule_set(imolecule)%lci)
305 : END IF
306 : END DO
307 12138 : DEALLOCATE (molecule_set)
308 :
309 : END IF
310 12138 : NULLIFY (molecule_set)
311 :
312 12138 : END SUBROUTINE deallocate_molecule_set
313 :
314 : ! **************************************************************************************************
315 : !> \brief Get components from a molecule data set.
316 : !> \param molecule ...
317 : !> \param molecule_kind ...
318 : !> \param lmi ...
319 : !> \param lci ...
320 : !> \param lg3x3 ...
321 : !> \param lg4x6 ...
322 : !> \param lcolv ...
323 : !> \param first_atom ...
324 : !> \param last_atom ...
325 : !> \param first_shell ...
326 : !> \param last_shell ...
327 : !> \date 29.08.2003
328 : !> \author Matthias Krack
329 : !> \version 1.0
330 : ! **************************************************************************************************
331 8735950 : SUBROUTINE get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, &
332 : first_atom, last_atom, first_shell, last_shell)
333 :
334 : TYPE(molecule_type), INTENT(IN) :: molecule
335 : TYPE(molecule_kind_type), OPTIONAL, POINTER :: molecule_kind
336 : TYPE(local_states_type), DIMENSION(:), OPTIONAL, &
337 : POINTER :: lmi
338 : TYPE(local_constraint_type), OPTIONAL, POINTER :: lci
339 : TYPE(local_g3x3_constraint_type), OPTIONAL, &
340 : POINTER :: lg3x3(:)
341 : TYPE(local_g4x6_constraint_type), OPTIONAL, &
342 : POINTER :: lg4x6(:)
343 : TYPE(local_colvar_constraint_type), DIMENSION(:), &
344 : OPTIONAL, POINTER :: lcolv
345 : INTEGER, OPTIONAL :: first_atom, last_atom, first_shell, &
346 : last_shell
347 :
348 8735950 : IF (PRESENT(first_atom)) first_atom = molecule%first_atom
349 8735950 : IF (PRESENT(last_atom)) last_atom = molecule%last_atom
350 8735950 : IF (PRESENT(first_shell)) first_shell = molecule%first_shell
351 8735950 : IF (PRESENT(last_shell)) last_shell = molecule%last_shell
352 8735950 : IF (PRESENT(molecule_kind)) molecule_kind => molecule%molecule_kind
353 8735950 : IF (PRESENT(lmi)) lmi => molecule%lmi
354 8735950 : IF (PRESENT(lci)) lci => molecule%lci
355 8735950 : IF (PRESENT(lcolv)) THEN
356 928472 : IF (ASSOCIATED(molecule%lci)) THEN
357 928472 : lcolv => molecule%lci%lcolv
358 : ELSE
359 0 : CPABORT("The pointer lci is not associated")
360 : END IF
361 : END IF
362 8735950 : IF (PRESENT(lg3x3)) THEN
363 1529202 : IF (ASSOCIATED(molecule%lci)) THEN
364 1529202 : lg3x3 => molecule%lci%lg3x3
365 : ELSE
366 0 : CPABORT("The pointer lci is not associated")
367 : END IF
368 : END IF
369 8735950 : IF (PRESENT(lg4x6)) THEN
370 885468 : IF (ASSOCIATED(molecule%lci)) THEN
371 885468 : lg4x6 => molecule%lci%lg4x6
372 : ELSE
373 0 : CPABORT("The pointer lci is not associated")
374 : END IF
375 : END IF
376 :
377 8735950 : END SUBROUTINE get_molecule
378 :
379 : ! **************************************************************************************************
380 : !> \brief Set a molecule data set.
381 : !> \param molecule ...
382 : !> \param molecule_kind ...
383 : !> \param lmi ...
384 : !> \param lci ...
385 : !> \param lcolv ...
386 : !> \param lg3x3 ...
387 : !> \param lg4x6 ...
388 : !> \date 29.08.2003
389 : !> \author Matthias Krack
390 : !> \version 1.0
391 : ! **************************************************************************************************
392 994390 : SUBROUTINE set_molecule(molecule, molecule_kind, lmi, lci, lcolv, lg3x3, lg4x6)
393 : TYPE(molecule_type), INTENT(INOUT) :: molecule
394 : TYPE(molecule_kind_type), OPTIONAL, POINTER :: molecule_kind
395 : TYPE(local_states_type), DIMENSION(:), OPTIONAL, &
396 : POINTER :: lmi
397 : TYPE(local_constraint_type), OPTIONAL, POINTER :: lci
398 : TYPE(local_colvar_constraint_type), DIMENSION(:), &
399 : OPTIONAL, POINTER :: lcolv
400 : TYPE(local_g3x3_constraint_type), OPTIONAL, &
401 : POINTER :: lg3x3(:)
402 : TYPE(local_g4x6_constraint_type), OPTIONAL, &
403 : POINTER :: lg4x6(:)
404 :
405 994390 : IF (PRESENT(molecule_kind)) molecule%molecule_kind => molecule_kind
406 994390 : IF (PRESENT(lmi)) molecule%lmi => lmi
407 994390 : IF (PRESENT(lci)) molecule%lci => lci
408 994390 : IF (PRESENT(lcolv)) THEN
409 2138 : IF (ASSOCIATED(molecule%lci)) THEN
410 2138 : molecule%lci%lcolv => lcolv
411 : ELSE
412 0 : CPABORT("The pointer lci is not associated")
413 : END IF
414 : END IF
415 994390 : IF (PRESENT(lg3x3)) THEN
416 36290 : IF (ASSOCIATED(molecule%lci)) THEN
417 36290 : molecule%lci%lg3x3 => lg3x3
418 : ELSE
419 0 : CPABORT("The pointer lci is not associated")
420 : END IF
421 : END IF
422 994390 : IF (PRESENT(lg4x6)) THEN
423 650 : IF (ASSOCIATED(molecule%lci)) THEN
424 650 : molecule%lci%lg4x6 => lg4x6
425 : ELSE
426 0 : CPABORT("The pointer lci is not associated")
427 : END IF
428 : END IF
429 :
430 994390 : END SUBROUTINE set_molecule
431 :
432 : ! **************************************************************************************************
433 : !> \brief Set a molecule data set.
434 : !> \param molecule_set ...
435 : !> \param first_atom ...
436 : !> \param last_atom ...
437 : !> \date 29.08.2003
438 : !> \author Matthias Krack
439 : !> \version 1.0
440 : ! **************************************************************************************************
441 12139 : SUBROUTINE set_molecule_set(molecule_set, first_atom, last_atom)
442 : TYPE(molecule_type), DIMENSION(:), INTENT(INOUT) :: molecule_set
443 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: first_atom, last_atom
444 :
445 : INTEGER :: imolecule
446 :
447 12139 : IF (PRESENT(first_atom)) THEN
448 12139 : IF (SIZE(first_atom) /= SIZE(molecule_set)) THEN
449 : CALL cp_abort(__LOCATION__, &
450 : "The sizes of first_atom and molecule_set "// &
451 0 : "are different")
452 : END IF
453 :
454 331825 : DO imolecule = 1, SIZE(molecule_set)
455 331825 : molecule_set(imolecule)%first_atom = first_atom(imolecule)
456 : END DO
457 : END IF
458 :
459 12139 : IF (PRESENT(last_atom)) THEN
460 12139 : IF (SIZE(last_atom) /= SIZE(molecule_set)) THEN
461 : CALL cp_abort(__LOCATION__, &
462 : "The sizes of last_atom and molecule_set "// &
463 0 : "are different")
464 : END IF
465 :
466 331825 : DO imolecule = 1, SIZE(molecule_set)
467 331825 : molecule_set(imolecule)%last_atom = last_atom(imolecule)
468 : END DO
469 : END IF
470 :
471 12139 : END SUBROUTINE set_molecule_set
472 :
473 : ! **************************************************************************************************
474 : !> \brief finds for each atom the molecule it belongs to
475 : !> \param molecule_set ...
476 : !> \param atom_to_mol ...
477 : ! **************************************************************************************************
478 484 : SUBROUTINE molecule_of_atom(molecule_set, atom_to_mol)
479 : TYPE(molecule_type), DIMENSION(:), INTENT(IN) :: molecule_set
480 : INTEGER, DIMENSION(:), INTENT(OUT) :: atom_to_mol
481 :
482 : INTEGER :: first_atom, iatom, imol, last_atom
483 :
484 3174 : DO imol = 1, SIZE(molecule_set)
485 2690 : CALL get_molecule(molecule=molecule_set(imol), first_atom=first_atom, last_atom=last_atom)
486 9276 : DO iatom = first_atom, last_atom
487 8792 : atom_to_mol(iatom) = imol
488 : END DO ! iatom
489 : END DO ! imol
490 :
491 484 : END SUBROUTINE molecule_of_atom
492 :
493 : ! **************************************************************************************************
494 : !> \brief returns information about molecules in the set.
495 : !> \param molecule_set ...
496 : !> \param atom_to_mol ...
497 : !> \param mol_to_first_atom ...
498 : !> \param mol_to_last_atom ...
499 : !> \param mol_to_nelectrons ...
500 : !> \param mol_to_nbasis ...
501 : !> \param mol_to_charge ...
502 : !> \param mol_to_multiplicity ...
503 : !> \par History
504 : !> 2011.06 created [Rustam Z Khaliullin]
505 : !> \author Rustam Z Khaliullin
506 : ! **************************************************************************************************
507 780 : SUBROUTINE get_molecule_set_info(molecule_set, atom_to_mol, mol_to_first_atom, &
508 780 : mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, &
509 260 : mol_to_multiplicity)
510 :
511 : TYPE(molecule_type), DIMENSION(:), INTENT(IN) :: molecule_set
512 : INTEGER, DIMENSION(:), INTENT(OUT), OPTIONAL :: atom_to_mol, mol_to_first_atom, &
513 : mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity
514 :
515 : INTEGER :: first_atom, iatom, imol, last_atom, &
516 : nbasis, nelec
517 : REAL(KIND=dp) :: charge
518 : TYPE(molecule_kind_type), POINTER :: imol_kind
519 :
520 1948 : DO imol = 1, SIZE(molecule_set)
521 :
522 : CALL get_molecule(molecule=molecule_set(imol), molecule_kind=imol_kind, &
523 1688 : first_atom=first_atom, last_atom=last_atom)
524 :
525 1688 : IF (PRESENT(mol_to_nelectrons)) THEN
526 822 : CALL get_molecule_kind(imol_kind, nelectron=nelec)
527 822 : mol_to_nelectrons(imol) = nelec
528 : END IF
529 :
530 1688 : IF (PRESENT(mol_to_multiplicity)) THEN
531 : ! RZK-warning: At the moment we can only get the total number
532 : ! of electrons (alpha+beta) and we do not have a way to get the multiplicity of mols.
533 : ! Therefore, the best we can do is to assume the singlet state for even number of electrons
534 : ! and doublet state for odd number of electrons (assume ne_alpha > ne_beta).
535 : ! The best way to implement a correct multiplicity subroutine in the future is to get
536 : ! the number of alpha and beta e- for each atom from init_atom_electronic_state. This way (as opposed to
537 : ! reading the multiplicities from file) the number of occupied and virtual orbitals
538 : ! will be consistent with atomic guess. A guess with broken symmetry will be easy to
539 : ! implement as well.
540 854 : CALL get_molecule_kind(imol_kind, nelectron=nelec)
541 854 : IF (MOD(nelec, 2) == 0) THEN
542 844 : mol_to_multiplicity(imol) = 1
543 : ELSE
544 10 : mol_to_multiplicity(imol) = 2
545 : END IF
546 : END IF
547 :
548 1688 : IF (PRESENT(mol_to_charge)) THEN
549 854 : CALL get_molecule_kind(imol_kind, charge=charge)
550 854 : mol_to_charge(imol) = NINT(charge)
551 : END IF
552 :
553 1688 : IF (PRESENT(mol_to_nbasis)) THEN
554 822 : CALL get_molecule_kind(imol_kind, nsgf=nbasis)
555 822 : mol_to_nbasis(imol) = nbasis
556 : END IF
557 :
558 1688 : IF (PRESENT(mol_to_first_atom)) THEN
559 1688 : mol_to_first_atom(imol) = first_atom
560 : END IF
561 :
562 1688 : IF (PRESENT(mol_to_last_atom)) THEN
563 1688 : mol_to_last_atom(imol) = last_atom
564 : END IF
565 :
566 1948 : IF (PRESENT(atom_to_mol)) THEN
567 2806 : DO iatom = first_atom, last_atom
568 2806 : atom_to_mol(iatom) = imol
569 : END DO ! iatom
570 : END IF
571 :
572 : END DO ! imol
573 :
574 260 : END SUBROUTINE get_molecule_set_info
575 :
576 : ! **************************************************************************************************
577 : !> \brief ...
578 : !> \param molecule_set ...
579 : !> \param atom_to_mol ...
580 : !> \param mol_to_first_atom ...
581 : !> \param mol_to_last_atom ...
582 : !> \param mol_to_nelectrons ...
583 : !> \param mol_to_nbasis ...
584 : !> \param mol_to_charge ...
585 : !> \param mol_to_multiplicity ...
586 : ! **************************************************************************************************
587 0 : SUBROUTINE get_domain_set_info(molecule_set, atom_to_mol, mol_to_first_atom, &
588 0 : mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, &
589 0 : mol_to_multiplicity)
590 :
591 : TYPE(molecule_type), DIMENSION(:), INTENT(IN) :: molecule_set
592 : INTEGER, DIMENSION(:), INTENT(OUT), OPTIONAL :: atom_to_mol, mol_to_first_atom, &
593 : mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity
594 :
595 : INTEGER :: first_atom, iatom, imol, last_atom, &
596 : nbasis, nelec
597 : REAL(KIND=dp) :: charge
598 : TYPE(molecule_kind_type), POINTER :: imol_kind
599 :
600 0 : DO imol = 1, SIZE(molecule_set)
601 :
602 : CALL get_molecule(molecule=molecule_set(imol), molecule_kind=imol_kind, &
603 0 : first_atom=first_atom, last_atom=last_atom)
604 :
605 0 : IF (PRESENT(mol_to_nelectrons)) THEN
606 0 : CALL get_molecule_kind(imol_kind, nelectron=nelec)
607 0 : mol_to_nelectrons(imol) = nelec
608 : END IF
609 :
610 0 : IF (PRESENT(mol_to_multiplicity)) THEN
611 : ! RZK-warning: At the moment we can only get the total number
612 : ! of electrons (alpha+beta) and we do not have a way to get the multiplicity of mols.
613 : ! Therefore, the best we can do is to assume the singlet state for even number of electrons
614 : ! and doublet state for odd number of electrons (assume ne_alpha > ne_beta).
615 : ! The best way to implement a correct multiplicity subroutine in the future is to get
616 : ! the number of alpha and beta e- for each atom from init_atom_electronic_state. This way (as opposed to
617 : ! reading the multiplicities from file) the number of occupied and virtual orbitals
618 : ! will be consistent with atomic guess. A guess with broken symmetry will be easy to
619 : ! implement as well.
620 0 : CALL get_molecule_kind(imol_kind, nelectron=nelec)
621 0 : IF (MOD(nelec, 2) == 0) THEN
622 0 : mol_to_multiplicity(imol) = 1
623 : ELSE
624 0 : mol_to_multiplicity(imol) = 2
625 : END IF
626 : END IF
627 :
628 0 : IF (PRESENT(mol_to_charge)) THEN
629 0 : CALL get_molecule_kind(imol_kind, charge=charge)
630 0 : mol_to_charge(imol) = NINT(charge)
631 : END IF
632 :
633 0 : IF (PRESENT(mol_to_nbasis)) THEN
634 0 : CALL get_molecule_kind(imol_kind, nsgf=nbasis)
635 0 : mol_to_nbasis(imol) = nbasis
636 : END IF
637 :
638 0 : IF (PRESENT(mol_to_first_atom)) THEN
639 0 : mol_to_first_atom(imol) = first_atom
640 : END IF
641 :
642 0 : IF (PRESENT(mol_to_last_atom)) THEN
643 0 : mol_to_last_atom(imol) = last_atom
644 : END IF
645 :
646 0 : IF (PRESENT(atom_to_mol)) THEN
647 0 : DO iatom = first_atom, last_atom
648 0 : atom_to_mol(iatom) = imol
649 : END DO ! iatom
650 : END IF
651 :
652 : END DO ! imol
653 :
654 0 : END SUBROUTINE get_domain_set_info
655 :
656 0 : END MODULE molecule_types
|
