(git:0761fd1)
Loading...
Searching...
No Matches
molecule_types.F
Go to the documentation of this file.
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! **************************************************************************************************
16MODULE 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
192CONTAINS
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 SUBROUTINE deallocate_global_constraint(gci)
203 TYPE(global_constraint_type), POINTER :: gci
204
205 INTEGER :: i
206
207 IF (ASSOCIATED(gci)) THEN
208 ! List of constraints
209 IF (ASSOCIATED(gci%colv_list)) THEN
210 DO i = 1, SIZE(gci%colv_list)
211 DEALLOCATE (gci%colv_list(i)%i_atoms)
212 END DO
213 DEALLOCATE (gci%colv_list)
214 END IF
215
216 IF (ASSOCIATED(gci%g3x3_list)) THEN
217 DEALLOCATE (gci%g3x3_list)
218 END IF
219
220 IF (ASSOCIATED(gci%g4x6_list)) THEN
221 DEALLOCATE (gci%g4x6_list)
222 END IF
223
224 ! Local information
225 IF (ASSOCIATED(gci%lcolv)) THEN
226 DO i = 1, SIZE(gci%lcolv)
227 CALL colvar_release(gci%lcolv(i)%colvar)
228 CALL colvar_release(gci%lcolv(i)%colvar_old)
229 END DO
230 DEALLOCATE (gci%lcolv)
231 END IF
232
233 IF (ASSOCIATED(gci%lg3x3)) THEN
234 DEALLOCATE (gci%lg3x3)
235 END IF
236
237 IF (ASSOCIATED(gci%lg4x6)) THEN
238 DEALLOCATE (gci%lg4x6)
239 END IF
240
241 IF (ASSOCIATED(gci%fixd_list)) THEN
242 DEALLOCATE (gci%fixd_list)
243 END IF
244
245 DEALLOCATE (gci)
246 END IF
247 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 SUBROUTINE allocate_molecule_set(molecule_set, nmolecule)
258 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
259 INTEGER, INTENT(IN) :: nmolecule
260
261 IF (ASSOCIATED(molecule_set)) CALL deallocate_molecule_set(molecule_set)
262
263 ALLOCATE (molecule_set(nmolecule))
264
265 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 SUBROUTINE deallocate_molecule_set(molecule_set)
275 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
276
277 INTEGER :: imolecule, j
278
279 IF (ASSOCIATED(molecule_set)) THEN
280
281 DO imolecule = 1, SIZE(molecule_set)
282 IF (ASSOCIATED(molecule_set(imolecule)%lmi)) THEN
283 DO j = 1, SIZE(molecule_set(imolecule)%lmi)
284 IF (ASSOCIATED(molecule_set(imolecule)%lmi(j)%states)) THEN
285 DEALLOCATE (molecule_set(imolecule)%lmi(j)%states)
286 END IF
287 END DO
288 DEALLOCATE (molecule_set(imolecule)%lmi)
289 END IF
290 IF (ASSOCIATED(molecule_set(imolecule)%lci)) THEN
291 IF (ASSOCIATED(molecule_set(imolecule)%lci%lcolv)) THEN
292 DO j = 1, SIZE(molecule_set(imolecule)%lci%lcolv)
293 CALL colvar_release(molecule_set(imolecule)%lci%lcolv(j)%colvar)
294 CALL colvar_release(molecule_set(imolecule)%lci%lcolv(j)%colvar_old)
295 END DO
296 DEALLOCATE (molecule_set(imolecule)%lci%lcolv)
297 END IF
298 IF (ASSOCIATED(molecule_set(imolecule)%lci%lg3x3)) THEN
299 DEALLOCATE (molecule_set(imolecule)%lci%lg3x3)
300 END IF
301 IF (ASSOCIATED(molecule_set(imolecule)%lci%lg4x6)) THEN
302 DEALLOCATE (molecule_set(imolecule)%lci%lg4x6)
303 END IF
304 DEALLOCATE (molecule_set(imolecule)%lci)
305 END IF
306 END DO
307 DEALLOCATE (molecule_set)
308
309 END IF
310 NULLIFY (molecule_set)
311
312 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 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 IF (PRESENT(first_atom)) first_atom = molecule%first_atom
349 IF (PRESENT(last_atom)) last_atom = molecule%last_atom
350 IF (PRESENT(first_shell)) first_shell = molecule%first_shell
351 IF (PRESENT(last_shell)) last_shell = molecule%last_shell
352 IF (PRESENT(molecule_kind)) molecule_kind => molecule%molecule_kind
353 IF (PRESENT(lmi)) lmi => molecule%lmi
354 IF (PRESENT(lci)) lci => molecule%lci
355 IF (PRESENT(lcolv)) THEN
356 IF (ASSOCIATED(molecule%lci)) THEN
357 lcolv => molecule%lci%lcolv
358 ELSE
359 cpabort("The pointer lci is not associated")
360 END IF
361 END IF
362 IF (PRESENT(lg3x3)) THEN
363 IF (ASSOCIATED(molecule%lci)) THEN
364 lg3x3 => molecule%lci%lg3x3
365 ELSE
366 cpabort("The pointer lci is not associated")
367 END IF
368 END IF
369 IF (PRESENT(lg4x6)) THEN
370 IF (ASSOCIATED(molecule%lci)) THEN
371 lg4x6 => molecule%lci%lg4x6
372 ELSE
373 cpabort("The pointer lci is not associated")
374 END IF
375 END IF
376
377 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 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 IF (PRESENT(molecule_kind)) molecule%molecule_kind => molecule_kind
406 IF (PRESENT(lmi)) molecule%lmi => lmi
407 IF (PRESENT(lci)) molecule%lci => lci
408 IF (PRESENT(lcolv)) THEN
409 IF (ASSOCIATED(molecule%lci)) THEN
410 molecule%lci%lcolv => lcolv
411 ELSE
412 cpabort("The pointer lci is not associated")
413 END IF
414 END IF
415 IF (PRESENT(lg3x3)) THEN
416 IF (ASSOCIATED(molecule%lci)) THEN
417 molecule%lci%lg3x3 => lg3x3
418 ELSE
419 cpabort("The pointer lci is not associated")
420 END IF
421 END IF
422 IF (PRESENT(lg4x6)) THEN
423 IF (ASSOCIATED(molecule%lci)) THEN
424 molecule%lci%lg4x6 => lg4x6
425 ELSE
426 cpabort("The pointer lci is not associated")
427 END IF
428 END IF
429
430 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 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 IF (PRESENT(first_atom)) THEN
448 IF (SIZE(first_atom) /= SIZE(molecule_set)) THEN
449 CALL cp_abort(__location__, &
450 "The sizes of first_atom and molecule_set "// &
451 "are different")
452 END IF
453
454 DO imolecule = 1, SIZE(molecule_set)
455 molecule_set(imolecule)%first_atom = first_atom(imolecule)
456 END DO
457 END IF
458
459 IF (PRESENT(last_atom)) THEN
460 IF (SIZE(last_atom) /= SIZE(molecule_set)) THEN
461 CALL cp_abort(__location__, &
462 "The sizes of last_atom and molecule_set "// &
463 "are different")
464 END IF
465
466 DO imolecule = 1, SIZE(molecule_set)
467 molecule_set(imolecule)%last_atom = last_atom(imolecule)
468 END DO
469 END IF
470
471 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 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 DO imol = 1, SIZE(molecule_set)
485 CALL get_molecule(molecule=molecule_set(imol), first_atom=first_atom, last_atom=last_atom)
486 DO iatom = first_atom, last_atom
487 atom_to_mol(iatom) = imol
488 END DO ! iatom
489 END DO ! imol
490
491 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 SUBROUTINE get_molecule_set_info(molecule_set, atom_to_mol, mol_to_first_atom, &
508 mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, &
509 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 DO imol = 1, SIZE(molecule_set)
521
522 CALL get_molecule(molecule=molecule_set(imol), molecule_kind=imol_kind, &
523 first_atom=first_atom, last_atom=last_atom)
524
525 IF (PRESENT(mol_to_nelectrons)) THEN
526 CALL get_molecule_kind(imol_kind, nelectron=nelec)
527 mol_to_nelectrons(imol) = nelec
528 END IF
529
530 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 CALL get_molecule_kind(imol_kind, nelectron=nelec)
541 IF (mod(nelec, 2) == 0) THEN
542 mol_to_multiplicity(imol) = 1
543 ELSE
544 mol_to_multiplicity(imol) = 2
545 END IF
546 END IF
547
548 IF (PRESENT(mol_to_charge)) THEN
549 CALL get_molecule_kind(imol_kind, charge=charge)
550 mol_to_charge(imol) = nint(charge)
551 END IF
552
553 IF (PRESENT(mol_to_nbasis)) THEN
554 CALL get_molecule_kind(imol_kind, nsgf=nbasis)
555 mol_to_nbasis(imol) = nbasis
556 END IF
557
558 IF (PRESENT(mol_to_first_atom)) THEN
559 mol_to_first_atom(imol) = first_atom
560 END IF
561
562 IF (PRESENT(mol_to_last_atom)) THEN
563 mol_to_last_atom(imol) = last_atom
564 END IF
565
566 IF (PRESENT(atom_to_mol)) THEN
567 DO iatom = first_atom, last_atom
568 atom_to_mol(iatom) = imol
569 END DO ! iatom
570 END IF
571
572 END DO ! imol
573
574 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 SUBROUTINE get_domain_set_info(molecule_set, atom_to_mol, mol_to_first_atom, &
588 mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, &
589 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 DO imol = 1, SIZE(molecule_set)
601
602 CALL get_molecule(molecule=molecule_set(imol), molecule_kind=imol_kind, &
603 first_atom=first_atom, last_atom=last_atom)
604
605 IF (PRESENT(mol_to_nelectrons)) THEN
606 CALL get_molecule_kind(imol_kind, nelectron=nelec)
607 mol_to_nelectrons(imol) = nelec
608 END IF
609
610 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 CALL get_molecule_kind(imol_kind, nelectron=nelec)
621 IF (mod(nelec, 2) == 0) THEN
622 mol_to_multiplicity(imol) = 1
623 ELSE
624 mol_to_multiplicity(imol) = 2
625 END IF
626 END IF
627
628 IF (PRESENT(mol_to_charge)) THEN
629 CALL get_molecule_kind(imol_kind, charge=charge)
630 mol_to_charge(imol) = nint(charge)
631 END IF
632
633 IF (PRESENT(mol_to_nbasis)) THEN
634 CALL get_molecule_kind(imol_kind, nsgf=nbasis)
635 mol_to_nbasis(imol) = nbasis
636 END IF
637
638 IF (PRESENT(mol_to_first_atom)) THEN
639 mol_to_first_atom(imol) = first_atom
640 END IF
641
642 IF (PRESENT(mol_to_last_atom)) THEN
643 mol_to_last_atom(imol) = last_atom
644 END IF
645
646 IF (PRESENT(atom_to_mol)) THEN
647 DO iatom = first_atom, last_atom
648 atom_to_mol(iatom) = imol
649 END DO ! iatom
650 END IF
651
652 END DO ! imol
653
654 END SUBROUTINE get_domain_set_info
655
656END MODULE molecule_types
colvar_types
Initialize the collective variables types.
Definition colvar_types.F:15
colvar_types::colvar_release
recursive subroutine, public colvar_release(colvar)
releases the memory that might have been allocated by the colvar
Definition colvar_types.F:1205
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:518
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::deallocate_molecule_set
subroutine, public deallocate_molecule_set(molecule_set)
Deallocate a molecule set.
Definition molecule_types.F:275
molecule_types::get_domain_set_info
subroutine, public get_domain_set_info(molecule_set, atom_to_mol, mol_to_first_atom, mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity)
...
Definition molecule_types.F:590
molecule_types::deallocate_global_constraint
subroutine, public deallocate_global_constraint(gci)
Deallocate a global constraint.
Definition molecule_types.F:203
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition molecule_types.F:333
molecule_types::allocate_molecule_set
subroutine, public allocate_molecule_set(molecule_set, nmolecule)
Allocate a molecule set.
Definition molecule_types.F:258
molecule_types::molecule_of_atom
subroutine, public molecule_of_atom(molecule_set, atom_to_mol)
finds for each atom the molecule it belongs to
Definition molecule_types.F:479
molecule_types::set_molecule_set
subroutine, public set_molecule_set(molecule_set, first_atom, last_atom)
Set a molecule data set.
Definition molecule_types.F:442
molecule_types::set_molecule
subroutine, public set_molecule(molecule, molecule_kind, lmi, lci, lcolv, lg3x3, lg4x6)
Set a molecule data set.
Definition molecule_types.F:393
molecule_types::get_molecule_set_info
subroutine, public get_molecule_set_info(molecule_set, atom_to_mol, mol_to_first_atom, mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity)
returns information about molecules in the set.
Definition molecule_types.F:510
colvar_types::colvar_counters
Definition colvar_types.F:398
colvar_types::colvar_type
parameters for a collective variable
Definition colvar_types.F:355
molecule_kind_types::colvar_constraint_type
Definition molecule_kind_types.F:112
molecule_kind_types::fixd_constraint_type
Definition molecule_kind_types.F:142
molecule_kind_types::g3x3_constraint_type
Definition molecule_kind_types.F:123
molecule_kind_types::g4x6_constraint_type
Definition molecule_kind_types.F:129
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:153
molecule_kind_types::vsite_constraint_type
Definition molecule_kind_types.F:136
molecule_types::global_constraint_type
Definition molecule_types.F:134
molecule_types::local_colvar_constraint_type
Definition molecule_types.F:40
molecule_types::local_constraint_type
Definition molecule_types.F:125
molecule_types::local_g3x3_constraint_type
Definition molecule_types.F:48
molecule_types::local_g4x6_constraint_type
Definition molecule_types.F:76
molecule_types::local_states_type
Definition molecule_types.F:120
molecule_types::molecule_type
Definition molecule_types.F:160