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 : INTEGER :: ntot = 0, &
137 : nrestraint = 0, &
138 : ng3x3 = 0, &
139 : ng3x3_restraint = 0, &
140 : ng4x6 = 0, &
141 : ng4x6_restraint = 0, &
142 : nvsite = 0, &
143 : nvsite_restraint = 0
144 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list => NULL()
145 : TYPE(colvar_constraint_type), DIMENSION(:), &
146 : POINTER :: colv_list => NULL()
147 : TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list => NULL()
148 : TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list => NULL()
149 : TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list => NULL()
150 : TYPE(local_colvar_constraint_type), &
151 : DIMENSION(:), POINTER :: lcolv => NULL()
152 : TYPE(local_g3x3_constraint_type), DIMENSION(:), &
153 : POINTER :: lg3x3 => NULL()
154 : TYPE(local_g4x6_constraint_type), DIMENSION(:), &
155 : POINTER :: lg4x6 => NULL()
156 : END TYPE global_constraint_type
157 :
158 : ! Molecule type
159 : TYPE molecule_type
160 : TYPE(molecule_kind_type), POINTER :: molecule_kind => NULL() ! pointer to molecule kind information
161 : TYPE(local_states_type), DIMENSION(:), POINTER :: lmi => NULL() ! local (spin)-states information
162 : TYPE(local_constraint_type), POINTER :: lci => NULL() ! local molecule constraint info
163 : INTEGER :: first_atom = 0 ! global index of first atom in molecule
164 : INTEGER :: last_atom = 0 ! global index of last atom in molecule
165 : INTEGER :: first_shell = 0 ! global index of first shell atom in molecule
166 : INTEGER :: last_shell = 0 ! global index of last shell atom in molecule
167 : END TYPE molecule_type
168 :
169 : ! Public data types
170 :
171 : PUBLIC :: local_colvar_constraint_type, &
172 : local_g3x3_constraint_type, &
173 : local_g4x6_constraint_type, &
174 : local_constraint_type, &
175 : local_states_type, &
176 : global_constraint_type, &
177 : molecule_type
178 :
179 : ! Public subroutines
180 :
181 : PUBLIC :: deallocate_global_constraint, &
182 : allocate_molecule_set, &
183 : deallocate_molecule_set, &
184 : get_molecule, &
185 : set_molecule, &
186 : set_molecule_set, &
187 : molecule_of_atom, &
188 : get_molecule_set_info, &
189 : get_domain_set_info
190 :
191 : CONTAINS
192 :
193 : ! **************************************************************************************************
194 : !> \brief Deallocate a global constraint.
195 : !> \param gci ...
196 : !> \par History
197 : !> 07.2003 created [fawzi]
198 : !> 01.2014 moved from cp_subsys_release() into separate routine.
199 : !> \author Ole Schuett
200 : ! **************************************************************************************************
201 11338 : SUBROUTINE deallocate_global_constraint(gci)
202 : TYPE(global_constraint_type), POINTER :: gci
203 :
204 : INTEGER :: i
205 :
206 11338 : IF (ASSOCIATED(gci)) THEN
207 : ! List of constraints
208 10788 : IF (ASSOCIATED(gci%colv_list)) THEN
209 110 : DO i = 1, SIZE(gci%colv_list)
210 110 : DEALLOCATE (gci%colv_list(i)%i_atoms)
211 : END DO
212 44 : DEALLOCATE (gci%colv_list)
213 : END IF
214 :
215 10788 : IF (ASSOCIATED(gci%g3x3_list)) &
216 4 : DEALLOCATE (gci%g3x3_list)
217 :
218 10788 : IF (ASSOCIATED(gci%g4x6_list)) &
219 4 : DEALLOCATE (gci%g4x6_list)
220 :
221 : ! Local information
222 10788 : IF (ASSOCIATED(gci%lcolv)) THEN
223 110 : DO i = 1, SIZE(gci%lcolv)
224 66 : CALL colvar_release(gci%lcolv(i)%colvar)
225 110 : CALL colvar_release(gci%lcolv(i)%colvar_old)
226 : END DO
227 44 : DEALLOCATE (gci%lcolv)
228 : END IF
229 :
230 10788 : IF (ASSOCIATED(gci%lg3x3)) &
231 4 : DEALLOCATE (gci%lg3x3)
232 :
233 10788 : IF (ASSOCIATED(gci%lg4x6)) &
234 4 : DEALLOCATE (gci%lg4x6)
235 :
236 10788 : IF (ASSOCIATED(gci%fixd_list)) &
237 2 : DEALLOCATE (gci%fixd_list)
238 :
239 10788 : DEALLOCATE (gci)
240 : END IF
241 11338 : END SUBROUTINE deallocate_global_constraint
242 :
243 : ! **************************************************************************************************
244 : !> \brief Allocate a molecule set.
245 : !> \param molecule_set ...
246 : !> \param nmolecule ...
247 : !> \date 29.08.2003
248 : !> \author Matthias Krack
249 : !> \version 1.0
250 : ! **************************************************************************************************
251 11338 : SUBROUTINE allocate_molecule_set(molecule_set, nmolecule)
252 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
253 : INTEGER, INTENT(IN) :: nmolecule
254 :
255 11338 : IF (ASSOCIATED(molecule_set)) CALL deallocate_molecule_set(molecule_set)
256 :
257 348958 : ALLOCATE (molecule_set(nmolecule))
258 :
259 11338 : END SUBROUTINE allocate_molecule_set
260 :
261 : ! **************************************************************************************************
262 : !> \brief Deallocate a molecule set.
263 : !> \param molecule_set ...
264 : !> \date 29.08.2003
265 : !> \author Matthias Krack
266 : !> \version 1.0
267 : ! **************************************************************************************************
268 11338 : SUBROUTINE deallocate_molecule_set(molecule_set)
269 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
270 :
271 : INTEGER :: imolecule, j
272 :
273 11338 : IF (ASSOCIATED(molecule_set)) THEN
274 :
275 326282 : DO imolecule = 1, SIZE(molecule_set)
276 314944 : IF (ASSOCIATED(molecule_set(imolecule)%lmi)) THEN
277 70 : DO j = 1, SIZE(molecule_set(imolecule)%lmi)
278 70 : IF (ASSOCIATED(molecule_set(imolecule)%lmi(j)%states)) THEN
279 40 : DEALLOCATE (molecule_set(imolecule)%lmi(j)%states)
280 : END IF
281 : END DO
282 30 : DEALLOCATE (molecule_set(imolecule)%lmi)
283 : END IF
284 326282 : IF (ASSOCIATED(molecule_set(imolecule)%lci)) THEN
285 43692 : IF (ASSOCIATED(molecule_set(imolecule)%lci%lcolv)) THEN
286 4336 : DO j = 1, SIZE(molecule_set(imolecule)%lci%lcolv)
287 2228 : CALL colvar_release(molecule_set(imolecule)%lci%lcolv(j)%colvar)
288 4336 : CALL colvar_release(molecule_set(imolecule)%lci%lcolv(j)%colvar_old)
289 : END DO
290 2108 : DEALLOCATE (molecule_set(imolecule)%lci%lcolv)
291 : END IF
292 43692 : IF (ASSOCIATED(molecule_set(imolecule)%lci%lg3x3)) THEN
293 36354 : DEALLOCATE (molecule_set(imolecule)%lci%lg3x3)
294 : END IF
295 43692 : IF (ASSOCIATED(molecule_set(imolecule)%lci%lg4x6)) THEN
296 650 : DEALLOCATE (molecule_set(imolecule)%lci%lg4x6)
297 : END IF
298 43692 : DEALLOCATE (molecule_set(imolecule)%lci)
299 : END IF
300 : END DO
301 11338 : DEALLOCATE (molecule_set)
302 :
303 : END IF
304 11338 : NULLIFY (molecule_set)
305 :
306 11338 : END SUBROUTINE deallocate_molecule_set
307 :
308 : ! **************************************************************************************************
309 : !> \brief Get components from a molecule data set.
310 : !> \param molecule ...
311 : !> \param molecule_kind ...
312 : !> \param lmi ...
313 : !> \param lci ...
314 : !> \param lg3x3 ...
315 : !> \param lg4x6 ...
316 : !> \param lcolv ...
317 : !> \param first_atom ...
318 : !> \param last_atom ...
319 : !> \param first_shell ...
320 : !> \param last_shell ...
321 : !> \date 29.08.2003
322 : !> \author Matthias Krack
323 : !> \version 1.0
324 : ! **************************************************************************************************
325 8688461 : SUBROUTINE get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, &
326 : first_atom, last_atom, first_shell, last_shell)
327 :
328 : TYPE(molecule_type), INTENT(IN) :: molecule
329 : TYPE(molecule_kind_type), OPTIONAL, POINTER :: molecule_kind
330 : TYPE(local_states_type), DIMENSION(:), OPTIONAL, &
331 : POINTER :: lmi
332 : TYPE(local_constraint_type), OPTIONAL, POINTER :: lci
333 : TYPE(local_g3x3_constraint_type), OPTIONAL, &
334 : POINTER :: lg3x3(:)
335 : TYPE(local_g4x6_constraint_type), OPTIONAL, &
336 : POINTER :: lg4x6(:)
337 : TYPE(local_colvar_constraint_type), DIMENSION(:), &
338 : OPTIONAL, POINTER :: lcolv
339 : INTEGER, OPTIONAL :: first_atom, last_atom, first_shell, &
340 : last_shell
341 :
342 8688461 : IF (PRESENT(first_atom)) first_atom = molecule%first_atom
343 8688461 : IF (PRESENT(last_atom)) last_atom = molecule%last_atom
344 8688461 : IF (PRESENT(first_shell)) first_shell = molecule%first_shell
345 8688461 : IF (PRESENT(last_shell)) last_shell = molecule%last_shell
346 8688461 : IF (PRESENT(molecule_kind)) molecule_kind => molecule%molecule_kind
347 8688461 : IF (PRESENT(lmi)) lmi => molecule%lmi
348 8688461 : IF (PRESENT(lci)) lci => molecule%lci
349 8688461 : IF (PRESENT(lcolv)) THEN
350 928471 : IF (ASSOCIATED(molecule%lci)) THEN
351 928471 : lcolv => molecule%lci%lcolv
352 : ELSE
353 0 : CPABORT("The pointer lci is not associated")
354 : END IF
355 : END IF
356 8688461 : IF (PRESENT(lg3x3)) THEN
357 1530904 : IF (ASSOCIATED(molecule%lci)) THEN
358 1530904 : lg3x3 => molecule%lci%lg3x3
359 : ELSE
360 0 : CPABORT("The pointer lci is not associated")
361 : END IF
362 : END IF
363 8688461 : IF (PRESENT(lg4x6)) THEN
364 885788 : IF (ASSOCIATED(molecule%lci)) THEN
365 885788 : lg4x6 => molecule%lci%lg4x6
366 : ELSE
367 0 : CPABORT("The pointer lci is not associated")
368 : END IF
369 : END IF
370 :
371 8688461 : END SUBROUTINE get_molecule
372 :
373 : ! **************************************************************************************************
374 : !> \brief Set a molecule data set.
375 : !> \param molecule ...
376 : !> \param molecule_kind ...
377 : !> \param lmi ...
378 : !> \param lci ...
379 : !> \param lcolv ...
380 : !> \param lg3x3 ...
381 : !> \param lg4x6 ...
382 : !> \date 29.08.2003
383 : !> \author Matthias Krack
384 : !> \version 1.0
385 : ! **************************************************************************************************
386 980204 : SUBROUTINE set_molecule(molecule, molecule_kind, lmi, lci, lcolv, lg3x3, lg4x6)
387 : TYPE(molecule_type), INTENT(INOUT) :: molecule
388 : TYPE(molecule_kind_type), OPTIONAL, POINTER :: molecule_kind
389 : TYPE(local_states_type), DIMENSION(:), OPTIONAL, &
390 : POINTER :: lmi
391 : TYPE(local_constraint_type), OPTIONAL, POINTER :: lci
392 : TYPE(local_colvar_constraint_type), DIMENSION(:), &
393 : OPTIONAL, POINTER :: lcolv
394 : TYPE(local_g3x3_constraint_type), OPTIONAL, &
395 : POINTER :: lg3x3(:)
396 : TYPE(local_g4x6_constraint_type), OPTIONAL, &
397 : POINTER :: lg4x6(:)
398 :
399 980204 : IF (PRESENT(molecule_kind)) molecule%molecule_kind => molecule_kind
400 980204 : IF (PRESENT(lmi)) molecule%lmi => lmi
401 980204 : IF (PRESENT(lci)) molecule%lci => lci
402 980204 : IF (PRESENT(lcolv)) THEN
403 2108 : IF (ASSOCIATED(molecule%lci)) THEN
404 2108 : molecule%lci%lcolv => lcolv
405 : ELSE
406 0 : CPABORT("The pointer lci is not associated")
407 : END IF
408 : END IF
409 980204 : IF (PRESENT(lg3x3)) THEN
410 36354 : IF (ASSOCIATED(molecule%lci)) THEN
411 36354 : molecule%lci%lg3x3 => lg3x3
412 : ELSE
413 0 : CPABORT("The pointer lci is not associated")
414 : END IF
415 : END IF
416 980204 : IF (PRESENT(lg4x6)) THEN
417 650 : IF (ASSOCIATED(molecule%lci)) THEN
418 650 : molecule%lci%lg4x6 => lg4x6
419 : ELSE
420 0 : CPABORT("The pointer lci is not associated")
421 : END IF
422 : END IF
423 :
424 980204 : END SUBROUTINE set_molecule
425 :
426 : ! **************************************************************************************************
427 : !> \brief Set a molecule data set.
428 : !> \param molecule_set ...
429 : !> \param first_atom ...
430 : !> \param last_atom ...
431 : !> \date 29.08.2003
432 : !> \author Matthias Krack
433 : !> \version 1.0
434 : ! **************************************************************************************************
435 11338 : SUBROUTINE set_molecule_set(molecule_set, first_atom, last_atom)
436 : TYPE(molecule_type), DIMENSION(:), INTENT(INOUT) :: molecule_set
437 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: first_atom, last_atom
438 :
439 : INTEGER :: imolecule
440 :
441 11338 : IF (PRESENT(first_atom)) THEN
442 11338 : IF (SIZE(first_atom) /= SIZE(molecule_set)) THEN
443 : CALL cp_abort(__LOCATION__, &
444 : "The sizes of first_atom and molecule_set "// &
445 0 : "are different")
446 : END IF
447 :
448 326282 : DO imolecule = 1, SIZE(molecule_set)
449 326282 : molecule_set(imolecule)%first_atom = first_atom(imolecule)
450 : END DO
451 : END IF
452 :
453 11338 : IF (PRESENT(last_atom)) THEN
454 11338 : IF (SIZE(last_atom) /= SIZE(molecule_set)) THEN
455 : CALL cp_abort(__LOCATION__, &
456 : "The sizes of last_atom and molecule_set "// &
457 0 : "are different")
458 : END IF
459 :
460 326282 : DO imolecule = 1, SIZE(molecule_set)
461 326282 : molecule_set(imolecule)%last_atom = last_atom(imolecule)
462 : END DO
463 : END IF
464 :
465 11338 : END SUBROUTINE set_molecule_set
466 :
467 : ! **************************************************************************************************
468 : !> \brief finds for each atom the molecule it belongs to
469 : !> \param molecule_set ...
470 : !> \param atom_to_mol ...
471 : ! **************************************************************************************************
472 536 : SUBROUTINE molecule_of_atom(molecule_set, atom_to_mol)
473 : TYPE(molecule_type), DIMENSION(:), INTENT(IN) :: molecule_set
474 : INTEGER, DIMENSION(:), INTENT(OUT) :: atom_to_mol
475 :
476 : INTEGER :: first_atom, iatom, imol, last_atom
477 :
478 3386 : DO imol = 1, SIZE(molecule_set)
479 2850 : CALL get_molecule(molecule=molecule_set(imol), first_atom=first_atom, last_atom=last_atom)
480 9728 : DO iatom = first_atom, last_atom
481 9192 : atom_to_mol(iatom) = imol
482 : END DO ! iatom
483 : END DO ! imol
484 :
485 536 : END SUBROUTINE molecule_of_atom
486 :
487 : ! **************************************************************************************************
488 : !> \brief returns information about molecules in the set.
489 : !> \param molecule_set ...
490 : !> \param atom_to_mol ...
491 : !> \param mol_to_first_atom ...
492 : !> \param mol_to_last_atom ...
493 : !> \param mol_to_nelectrons ...
494 : !> \param mol_to_nbasis ...
495 : !> \param mol_to_charge ...
496 : !> \param mol_to_multiplicity ...
497 : !> \par History
498 : !> 2011.06 created [Rustam Z Khaliullin]
499 : !> \author Rustam Z Khaliullin
500 : ! **************************************************************************************************
501 780 : SUBROUTINE get_molecule_set_info(molecule_set, atom_to_mol, mol_to_first_atom, &
502 780 : mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, &
503 260 : mol_to_multiplicity)
504 :
505 : TYPE(molecule_type), DIMENSION(:), INTENT(IN) :: molecule_set
506 : INTEGER, DIMENSION(:), INTENT(OUT), OPTIONAL :: atom_to_mol, mol_to_first_atom, &
507 : mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity
508 :
509 : INTEGER :: first_atom, iatom, imol, last_atom, &
510 : nbasis, nelec
511 : REAL(KIND=dp) :: charge
512 : TYPE(molecule_kind_type), POINTER :: imol_kind
513 :
514 1948 : DO imol = 1, SIZE(molecule_set)
515 :
516 : CALL get_molecule(molecule=molecule_set(imol), molecule_kind=imol_kind, &
517 1688 : first_atom=first_atom, last_atom=last_atom)
518 :
519 1688 : IF (PRESENT(mol_to_nelectrons)) THEN
520 822 : CALL get_molecule_kind(imol_kind, nelectron=nelec)
521 822 : mol_to_nelectrons(imol) = nelec
522 : END IF
523 :
524 1688 : IF (PRESENT(mol_to_multiplicity)) THEN
525 : ! RZK-warning: At the moment we can only get the total number
526 : ! of electrons (alpha+beta) and we do not have a way to get the multiplicity of mols.
527 : ! Therefore, the best we can do is to assume the singlet state for even number of electrons
528 : ! and doublet state for odd number of electrons (assume ne_alpha > ne_beta).
529 : ! The best way to implement a correct multiplicity subroutine in the future is to get
530 : ! the number of alpha and beta e- for each atom from init_atom_electronic_state. This way (as opposed to
531 : ! reading the multiplicities from file) the number of occupied and virtual orbitals
532 : ! will be consistent with atomic guess. A guess with broken symmetry will be easy to
533 : ! implement as well.
534 854 : CALL get_molecule_kind(imol_kind, nelectron=nelec)
535 854 : IF (MOD(nelec, 2) == 0) THEN
536 844 : mol_to_multiplicity(imol) = 1
537 : ELSE
538 10 : mol_to_multiplicity(imol) = 2
539 : END IF
540 : END IF
541 :
542 1688 : IF (PRESENT(mol_to_charge)) THEN
543 854 : CALL get_molecule_kind(imol_kind, charge=charge)
544 854 : mol_to_charge(imol) = NINT(charge)
545 : END IF
546 :
547 1688 : IF (PRESENT(mol_to_nbasis)) THEN
548 822 : CALL get_molecule_kind(imol_kind, nsgf=nbasis)
549 822 : mol_to_nbasis(imol) = nbasis
550 : END IF
551 :
552 1688 : IF (PRESENT(mol_to_first_atom)) THEN
553 1688 : mol_to_first_atom(imol) = first_atom
554 : END IF
555 :
556 1688 : IF (PRESENT(mol_to_last_atom)) THEN
557 1688 : mol_to_last_atom(imol) = last_atom
558 : END IF
559 :
560 1948 : IF (PRESENT(atom_to_mol)) THEN
561 2806 : DO iatom = first_atom, last_atom
562 2806 : atom_to_mol(iatom) = imol
563 : END DO ! iatom
564 : END IF
565 :
566 : END DO ! imol
567 :
568 260 : END SUBROUTINE get_molecule_set_info
569 :
570 : ! **************************************************************************************************
571 : !> \brief ...
572 : !> \param molecule_set ...
573 : !> \param atom_to_mol ...
574 : !> \param mol_to_first_atom ...
575 : !> \param mol_to_last_atom ...
576 : !> \param mol_to_nelectrons ...
577 : !> \param mol_to_nbasis ...
578 : !> \param mol_to_charge ...
579 : !> \param mol_to_multiplicity ...
580 : ! **************************************************************************************************
581 0 : SUBROUTINE get_domain_set_info(molecule_set, atom_to_mol, mol_to_first_atom, &
582 0 : mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, &
583 0 : mol_to_multiplicity)
584 :
585 : TYPE(molecule_type), DIMENSION(:), INTENT(IN) :: molecule_set
586 : INTEGER, DIMENSION(:), INTENT(OUT), OPTIONAL :: atom_to_mol, mol_to_first_atom, &
587 : mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity
588 :
589 : INTEGER :: first_atom, iatom, imol, last_atom, &
590 : nbasis, nelec
591 : REAL(KIND=dp) :: charge
592 : TYPE(molecule_kind_type), POINTER :: imol_kind
593 :
594 0 : DO imol = 1, SIZE(molecule_set)
595 :
596 : CALL get_molecule(molecule=molecule_set(imol), molecule_kind=imol_kind, &
597 0 : first_atom=first_atom, last_atom=last_atom)
598 :
599 0 : IF (PRESENT(mol_to_nelectrons)) THEN
600 0 : CALL get_molecule_kind(imol_kind, nelectron=nelec)
601 0 : mol_to_nelectrons(imol) = nelec
602 : END IF
603 :
604 0 : IF (PRESENT(mol_to_multiplicity)) THEN
605 : ! RZK-warning: At the moment we can only get the total number
606 : ! of electrons (alpha+beta) and we do not have a way to get the multiplicity of mols.
607 : ! Therefore, the best we can do is to assume the singlet state for even number of electrons
608 : ! and doublet state for odd number of electrons (assume ne_alpha > ne_beta).
609 : ! The best way to implement a correct multiplicity subroutine in the future is to get
610 : ! the number of alpha and beta e- for each atom from init_atom_electronic_state. This way (as opposed to
611 : ! reading the multiplicities from file) the number of occupied and virtual orbitals
612 : ! will be consistent with atomic guess. A guess with broken symmetry will be easy to
613 : ! implement as well.
614 0 : CALL get_molecule_kind(imol_kind, nelectron=nelec)
615 0 : IF (MOD(nelec, 2) == 0) THEN
616 0 : mol_to_multiplicity(imol) = 1
617 : ELSE
618 0 : mol_to_multiplicity(imol) = 2
619 : END IF
620 : END IF
621 :
622 0 : IF (PRESENT(mol_to_charge)) THEN
623 0 : CALL get_molecule_kind(imol_kind, charge=charge)
624 0 : mol_to_charge(imol) = NINT(charge)
625 : END IF
626 :
627 0 : IF (PRESENT(mol_to_nbasis)) THEN
628 0 : CALL get_molecule_kind(imol_kind, nsgf=nbasis)
629 0 : mol_to_nbasis(imol) = nbasis
630 : END IF
631 :
632 0 : IF (PRESENT(mol_to_first_atom)) THEN
633 0 : mol_to_first_atom(imol) = first_atom
634 : END IF
635 :
636 0 : IF (PRESENT(mol_to_last_atom)) THEN
637 0 : mol_to_last_atom(imol) = last_atom
638 : END IF
639 :
640 0 : IF (PRESENT(atom_to_mol)) THEN
641 0 : DO iatom = first_atom, last_atom
642 0 : atom_to_mol(iatom) = imol
643 : END DO ! iatom
644 : END IF
645 :
646 : END DO ! imol
647 :
648 0 : END SUBROUTINE get_domain_set_info
649 :
650 0 : END MODULE molecule_types
|
