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