(git:374b731)
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-2024 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 MK (29.08.2003)
15! **************************************************************************************************
17
18 USE colvar_types, ONLY: colvar_counters,&
21 USE kinds, ONLY: dp
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! *** Data types ***
40! **************************************************************************************************
42 LOGICAL :: init = .false.
43 TYPE(colvar_type), POINTER :: colvar => null()
44 TYPE(colvar_type), POINTER :: colvar_old => null()
45 REAL(kind=dp) :: lambda = 0.0_dp, sigma = 0.0_dp
47
48! **************************************************************************************************
50 LOGICAL :: init = .false.
51 REAL(kind=dp) :: scale = 0.0_dp, scale_old = 0.0_dp, &
52 imass1 = 0.0_dp, imass2 = 0.0_dp, imass3 = 0.0_dp
53 REAL(kind=dp), DIMENSION(3) :: fa = 0.0_dp, fb = 0.0_dp, fc = 0.0_dp, &
54 f_roll1 = 0.0_dp, f_roll2 = 0.0_dp, f_roll3 = 0.0_dp, &
55 ra_old = 0.0_dp, rb_old = 0.0_dp, rc_old = 0.0_dp, &
56 va = 0.0_dp, vb = 0.0_dp, vc = 0.0_dp, &
57 lambda = 0.0_dp, del_lambda = 0.0_dp, lambda_old = 0.0_dp, &
58 r0_12 = 0.0_dp, r0_13 = 0.0_dp, r0_23 = 0.0_dp
59 REAL(kind=dp), DIMENSION(3, 3) :: amat = 0.0_dp
61
62! **************************************************************************************************
64 LOGICAL :: init = .false.
65 REAL(kind=dp) :: scale = 0.0_dp, scale_old = 0.0_dp, imass1 = 0.0_dp, &
66 imass2 = 0.0_dp, imass3 = 0.0_dp, imass4 = 0.0_dp
67 REAL(kind=dp), DIMENSION(3) :: fa = 0.0_dp, fb = 0.0_dp, fc = 0.0_dp, fd = 0.0_dp, fe = 0.0_dp, ff = 0.0_dp, &
68 f_roll1 = 0.0_dp, f_roll2 = 0.0_dp, f_roll3 = 0.0_dp, &
69 f_roll4 = 0.0_dp, f_roll5 = 0.0_dp, f_roll6 = 0.0_dp, &
70 ra_old = 0.0_dp, rb_old = 0.0_dp, rc_old = 0.0_dp, &
71 rd_old = 0.0_dp, re_old = 0.0_dp, rf_old = 0.0_dp, &
72 va = 0.0_dp, vb = 0.0_dp, vc = 0.0_dp, vd = 0.0_dp, ve = 0.0_dp, vf = 0.0_dp, &
73 r0_12 = 0.0_dp, r0_13 = 0.0_dp, r0_14 = 0.0_dp, &
74 r0_23 = 0.0_dp, r0_24 = 0.0_dp, r0_34 = 0.0_dp
75 REAL(kind=dp), DIMENSION(6) :: lambda = 0.0_dp, del_lambda = 0.0_dp, lambda_old = 0.0_dp
76 REAL(kind=dp), DIMENSION(6, 6) :: amat = 0.0_dp
78
79! **************************************************************************************************
81 INTEGER, DIMENSION(:), POINTER :: states => null() ! indices of Kohn-Sham states for molecule
82 INTEGER :: nstates = 0 ! Kohn-Sham states for molecule
83 END TYPE local_states_type
84
85! **************************************************************************************************
87 TYPE(local_colvar_constraint_type), DIMENSION(:), POINTER :: lcolv => null()
88 TYPE(local_g3x3_constraint_type), DIMENSION(:), POINTER :: lg3x3 => null()
89 TYPE(local_g4x6_constraint_type), DIMENSION(:), POINTER :: lg4x6 => null()
91
92! **************************************************************************************************
95 INTEGER :: ntot = 0, nrestraint = 0
96 INTEGER :: ng3x3 = 0, ng3x3_restraint = 0
97 INTEGER :: ng4x6 = 0, ng4x6_restraint = 0
98 INTEGER :: nvsite = 0, nvsite_restraint = 0
99 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list => null()
100 TYPE(colvar_constraint_type), DIMENSION(:), POINTER :: colv_list => null()
101 TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list => null()
102 TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list => null()
103 TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list => null()
104 TYPE(local_colvar_constraint_type), DIMENSION(:), POINTER :: lcolv => null()
105 TYPE(local_g3x3_constraint_type), DIMENSION(:), POINTER :: lg3x3 => null()
106 TYPE(local_g4x6_constraint_type), DIMENSION(:), POINTER :: lg4x6 => null()
108
109! **************************************************************************************************
111 TYPE(molecule_kind_type), POINTER :: molecule_kind => null() ! pointer to molecule kind information
112 TYPE(local_states_type), DIMENSION(:), POINTER :: lmi => null() ! local (spin)-states information
113 TYPE(local_constraint_type), POINTER :: lci => null() ! local molecule constraint info
114 INTEGER :: first_atom = 0 ! global index of first atom in molecule
115 INTEGER :: last_atom = 0 ! global index of last atom in molecule
116 INTEGER :: first_shell = 0 ! global index of first shell atom in molecule
117 INTEGER :: last_shell = 0 ! global index of last shell atom in molecule
118 END TYPE molecule_type
119
120! *** Public data types ***
121
129
130! *** Public subroutines ***
131
135 get_molecule, &
136 set_molecule, &
140
141CONTAINS
142
143! **************************************************************************************************
144!> \brief Deallocate a global constraint.
145!> \param gci ...
146!> \par History
147!> 07.2003 created [fawzi]
148!> 01.2014 moved from cp_subsys_release() into separate routine.
149!> \author Ole Schuett
150! **************************************************************************************************
152 TYPE(global_constraint_type), POINTER :: gci
153
154 INTEGER :: i
155
156 IF (ASSOCIATED(gci)) THEN
157 ! List of constraints
158 IF (ASSOCIATED(gci%colv_list)) THEN
159 DO i = 1, SIZE(gci%colv_list)
160 DEALLOCATE (gci%colv_list(i)%i_atoms)
161 END DO
162 DEALLOCATE (gci%colv_list)
163 END IF
164
165 IF (ASSOCIATED(gci%g3x3_list)) &
166 DEALLOCATE (gci%g3x3_list)
167
168 IF (ASSOCIATED(gci%g4x6_list)) &
169 DEALLOCATE (gci%g4x6_list)
170
171 ! Local information
172 IF (ASSOCIATED(gci%lcolv)) THEN
173 DO i = 1, SIZE(gci%lcolv)
174 CALL colvar_release(gci%lcolv(i)%colvar)
175 CALL colvar_release(gci%lcolv(i)%colvar_old)
176 END DO
177 DEALLOCATE (gci%lcolv)
178 END IF
179
180 IF (ASSOCIATED(gci%lg3x3)) &
181 DEALLOCATE (gci%lg3x3)
182
183 IF (ASSOCIATED(gci%lg4x6)) &
184 DEALLOCATE (gci%lg4x6)
185
186 IF (ASSOCIATED(gci%fixd_list)) &
187 DEALLOCATE (gci%fixd_list)
188
189 DEALLOCATE (gci)
190 END IF
191 END SUBROUTINE deallocate_global_constraint
192
193! **************************************************************************************************
194!> \brief Allocate a molecule set.
195!> \param molecule_set ...
196!> \param nmolecule ...
197!> \date 29.08.2003
198!> \author MK
199!> \version 1.0
200! **************************************************************************************************
201 SUBROUTINE allocate_molecule_set(molecule_set, nmolecule)
202 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
203 INTEGER, INTENT(IN) :: nmolecule
204
205 IF (ASSOCIATED(molecule_set)) CALL deallocate_molecule_set(molecule_set)
206
207 ALLOCATE (molecule_set(nmolecule))
208
209 END SUBROUTINE allocate_molecule_set
210
211! **************************************************************************************************
212!> \brief Deallocate a molecule set.
213!> \param molecule_set ...
214!> \date 29.08.2003
215!> \author MK
216!> \version 1.0
217! **************************************************************************************************
218 SUBROUTINE deallocate_molecule_set(molecule_set)
219 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
220
221 INTEGER :: imolecule, j
222
223 IF (ASSOCIATED(molecule_set)) THEN
224
225 DO imolecule = 1, SIZE(molecule_set)
226 IF (ASSOCIATED(molecule_set(imolecule)%lmi)) THEN
227 DO j = 1, SIZE(molecule_set(imolecule)%lmi)
228 IF (ASSOCIATED(molecule_set(imolecule)%lmi(j)%states)) THEN
229 DEALLOCATE (molecule_set(imolecule)%lmi(j)%states)
230 END IF
231 END DO
232 DEALLOCATE (molecule_set(imolecule)%lmi)
233 END IF
234 IF (ASSOCIATED(molecule_set(imolecule)%lci)) THEN
235 IF (ASSOCIATED(molecule_set(imolecule)%lci%lcolv)) THEN
236 DO j = 1, SIZE(molecule_set(imolecule)%lci%lcolv)
237 CALL colvar_release(molecule_set(imolecule)%lci%lcolv(j)%colvar)
238 CALL colvar_release(molecule_set(imolecule)%lci%lcolv(j)%colvar_old)
239 END DO
240 DEALLOCATE (molecule_set(imolecule)%lci%lcolv)
241 END IF
242 IF (ASSOCIATED(molecule_set(imolecule)%lci%lg3x3)) THEN
243 DEALLOCATE (molecule_set(imolecule)%lci%lg3x3)
244 END IF
245 IF (ASSOCIATED(molecule_set(imolecule)%lci%lg4x6)) THEN
246 DEALLOCATE (molecule_set(imolecule)%lci%lg4x6)
247 END IF
248 DEALLOCATE (molecule_set(imolecule)%lci)
249 END IF
250 END DO
251 DEALLOCATE (molecule_set)
252
253 END IF
254 NULLIFY (molecule_set)
255
256 END SUBROUTINE deallocate_molecule_set
257
258! **************************************************************************************************
259!> \brief Get components from a molecule data set.
260!> \param molecule ...
261!> \param molecule_kind ...
262!> \param lmi ...
263!> \param lci ...
264!> \param lg3x3 ...
265!> \param lg4x6 ...
266!> \param lcolv ...
267!> \param first_atom ...
268!> \param last_atom ...
269!> \param first_shell ...
270!> \param last_shell ...
271!> \date 29.08.2003
272!> \author MK
273!> \version 1.0
274! **************************************************************************************************
275 SUBROUTINE get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, &
276 first_atom, last_atom, first_shell, last_shell)
277
278 TYPE(molecule_type), INTENT(IN) :: molecule
279 TYPE(molecule_kind_type), OPTIONAL, POINTER :: molecule_kind
280 TYPE(local_states_type), DIMENSION(:), OPTIONAL, &
281 POINTER :: lmi
282 TYPE(local_constraint_type), OPTIONAL, POINTER :: lci
283 TYPE(local_g3x3_constraint_type), OPTIONAL, &
284 POINTER :: lg3x3(:)
285 TYPE(local_g4x6_constraint_type), OPTIONAL, &
286 POINTER :: lg4x6(:)
287 TYPE(local_colvar_constraint_type), DIMENSION(:), &
288 OPTIONAL, POINTER :: lcolv
289 INTEGER, OPTIONAL :: first_atom, last_atom, first_shell, &
290 last_shell
291
292 IF (PRESENT(first_atom)) first_atom = molecule%first_atom
293 IF (PRESENT(last_atom)) last_atom = molecule%last_atom
294 IF (PRESENT(first_shell)) first_shell = molecule%first_shell
295 IF (PRESENT(last_shell)) last_shell = molecule%last_shell
296 IF (PRESENT(molecule_kind)) molecule_kind => molecule%molecule_kind
297 IF (PRESENT(lmi)) lmi => molecule%lmi
298 IF (PRESENT(lci)) lci => molecule%lci
299 IF (PRESENT(lcolv)) THEN
300 IF (ASSOCIATED(molecule%lci)) THEN
301 lcolv => molecule%lci%lcolv
302 ELSE
303 cpabort("The pointer lci is not associated")
304 END IF
305 END IF
306 IF (PRESENT(lg3x3)) THEN
307 IF (ASSOCIATED(molecule%lci)) THEN
308 lg3x3 => molecule%lci%lg3x3
309 ELSE
310 cpabort("The pointer lci is not associated")
311 END IF
312 END IF
313 IF (PRESENT(lg4x6)) THEN
314 IF (ASSOCIATED(molecule%lci)) THEN
315 lg4x6 => molecule%lci%lg4x6
316 ELSE
317 cpabort("The pointer lci is not associated")
318 END IF
319 END IF
320
321 END SUBROUTINE get_molecule
322
323! **************************************************************************************************
324!> \brief Set a molecule data set.
325!> \param molecule ...
326!> \param molecule_kind ...
327!> \param lmi ...
328!> \param lci ...
329!> \param lcolv ...
330!> \param lg3x3 ...
331!> \param lg4x6 ...
332!> \date 29.08.2003
333!> \author MK
334!> \version 1.0
335! **************************************************************************************************
336 SUBROUTINE set_molecule(molecule, molecule_kind, lmi, lci, lcolv, lg3x3, lg4x6)
337 TYPE(molecule_type), INTENT(INOUT) :: molecule
338 TYPE(molecule_kind_type), OPTIONAL, POINTER :: molecule_kind
339 TYPE(local_states_type), DIMENSION(:), OPTIONAL, &
340 POINTER :: lmi
341 TYPE(local_constraint_type), OPTIONAL, POINTER :: lci
342 TYPE(local_colvar_constraint_type), DIMENSION(:), &
343 OPTIONAL, POINTER :: lcolv
344 TYPE(local_g3x3_constraint_type), OPTIONAL, &
345 POINTER :: lg3x3(:)
346 TYPE(local_g4x6_constraint_type), OPTIONAL, &
347 POINTER :: lg4x6(:)
348
349 IF (PRESENT(molecule_kind)) molecule%molecule_kind => molecule_kind
350 IF (PRESENT(lmi)) molecule%lmi => lmi
351 IF (PRESENT(lci)) molecule%lci => lci
352 IF (PRESENT(lcolv)) THEN
353 IF (ASSOCIATED(molecule%lci)) THEN
354 molecule%lci%lcolv => lcolv
355 ELSE
356 cpabort("The pointer lci is not associated")
357 END IF
358 END IF
359 IF (PRESENT(lg3x3)) THEN
360 IF (ASSOCIATED(molecule%lci)) THEN
361 molecule%lci%lg3x3 => lg3x3
362 ELSE
363 cpabort("The pointer lci is not associated")
364 END IF
365 END IF
366 IF (PRESENT(lg4x6)) THEN
367 IF (ASSOCIATED(molecule%lci)) THEN
368 molecule%lci%lg4x6 => lg4x6
369 ELSE
370 cpabort("The pointer lci is not associated")
371 END IF
372 END IF
373
374 END SUBROUTINE set_molecule
375
376! **************************************************************************************************
377!> \brief Set a molecule data set.
378!> \param molecule_set ...
379!> \param first_atom ...
380!> \param last_atom ...
381!> \date 29.08.2003
382!> \author MK
383!> \version 1.0
384! **************************************************************************************************
385 SUBROUTINE set_molecule_set(molecule_set, first_atom, last_atom)
386 TYPE(molecule_type), DIMENSION(:), INTENT(INOUT) :: molecule_set
387 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: first_atom, last_atom
388
389 INTEGER :: imolecule
390
391 IF (PRESENT(first_atom)) THEN
392 IF (SIZE(first_atom) /= SIZE(molecule_set)) THEN
393 CALL cp_abort(__location__, &
394 "The sizes of first_atom and molecule_set "// &
395 "are different")
396 END IF
397
398 DO imolecule = 1, SIZE(molecule_set)
399 molecule_set(imolecule)%first_atom = first_atom(imolecule)
400 END DO
401 END IF
402
403 IF (PRESENT(last_atom)) THEN
404 IF (SIZE(last_atom) /= SIZE(molecule_set)) THEN
405 CALL cp_abort(__location__, &
406 "The sizes of last_atom and molecule_set "// &
407 "are different")
408 END IF
409
410 DO imolecule = 1, SIZE(molecule_set)
411 molecule_set(imolecule)%last_atom = last_atom(imolecule)
412 END DO
413 END IF
414
415 END SUBROUTINE set_molecule_set
416
417! **************************************************************************************************
418!> \brief finds for each atom the molecule it belongs to
419!> \param molecule_set ...
420!> \param atom_to_mol ...
421! **************************************************************************************************
422 SUBROUTINE molecule_of_atom(molecule_set, atom_to_mol)
423 TYPE(molecule_type), DIMENSION(:), INTENT(IN) :: molecule_set
424 INTEGER, DIMENSION(:), INTENT(OUT) :: atom_to_mol
425
426 INTEGER :: first_atom, iatom, imol, last_atom
427
428 DO imol = 1, SIZE(molecule_set)
429 CALL get_molecule(molecule=molecule_set(imol), first_atom=first_atom, last_atom=last_atom)
430 DO iatom = first_atom, last_atom
431 atom_to_mol(iatom) = imol
432 END DO ! iatom
433 END DO ! imol
434
435 END SUBROUTINE molecule_of_atom
436
437! **************************************************************************************************
438!> \brief returns information about molecules in the set.
439!> \param molecule_set ...
440!> \param atom_to_mol ...
441!> \param mol_to_first_atom ...
442!> \param mol_to_last_atom ...
443!> \param mol_to_nelectrons ...
444!> \param mol_to_nbasis ...
445!> \param mol_to_charge ...
446!> \param mol_to_multiplicity ...
447!> \par History
448!> 2011.06 created [Rustam Z Khaliullin]
449!> \author Rustam Z Khaliullin
450! **************************************************************************************************
451 SUBROUTINE get_molecule_set_info(molecule_set, atom_to_mol, mol_to_first_atom, &
452 mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, &
453 mol_to_multiplicity)
454
455 TYPE(molecule_type), DIMENSION(:), INTENT(IN) :: molecule_set
456 INTEGER, DIMENSION(:), INTENT(OUT), OPTIONAL :: atom_to_mol, mol_to_first_atom, &
457 mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity
458
459 INTEGER :: first_atom, iatom, imol, last_atom, &
460 nbasis, nelec
461 REAL(kind=dp) :: charge
462 TYPE(molecule_kind_type), POINTER :: imol_kind
463
464 DO imol = 1, SIZE(molecule_set)
465
466 CALL get_molecule(molecule=molecule_set(imol), molecule_kind=imol_kind, &
467 first_atom=first_atom, last_atom=last_atom)
468
469 IF (PRESENT(mol_to_nelectrons)) THEN
470 CALL get_molecule_kind(imol_kind, nelectron=nelec)
471 mol_to_nelectrons(imol) = nelec
472 END IF
473
474 IF (PRESENT(mol_to_multiplicity)) THEN
475 ! RZK-warning: At the moment we can only get the total number
476 ! of electrons (alpha+beta) and we do not have a way to get the multiplicity of mols.
477 ! Therefore, the best we can do is to assume the singlet state for even number of electrons
478 ! and doublet state for odd number of electrons (assume ne_alpha > ne_beta).
479 ! The best way to implement a correct multiplicity subroutine in the future is to get
480 ! the number of alpha and beta e- for each atom from init_atom_electronic_state. This way (as opposed to
481 ! reading the multiplicities from file) the number of occupied and virtual orbitals
482 ! will be consistent with atomic guess. A guess with broken symmetry will be easy to
483 ! implement as well.
484 CALL get_molecule_kind(imol_kind, nelectron=nelec)
485 IF (mod(nelec, 2) .EQ. 0) THEN
486 mol_to_multiplicity(imol) = 1
487 ELSE
488 mol_to_multiplicity(imol) = 2
489 END IF
490 END IF
491
492 IF (PRESENT(mol_to_charge)) THEN
493 CALL get_molecule_kind(imol_kind, charge=charge)
494 mol_to_charge(imol) = nint(charge)
495 END IF
496
497 IF (PRESENT(mol_to_nbasis)) THEN
498 CALL get_molecule_kind(imol_kind, nsgf=nbasis)
499 mol_to_nbasis(imol) = nbasis
500 END IF
501
502 IF (PRESENT(mol_to_first_atom)) THEN
503 mol_to_first_atom(imol) = first_atom
504 END IF
505
506 IF (PRESENT(mol_to_last_atom)) THEN
507 mol_to_last_atom(imol) = last_atom
508 END IF
509
510 IF (PRESENT(atom_to_mol)) THEN
511 DO iatom = first_atom, last_atom
512 atom_to_mol(iatom) = imol
513 END DO ! iatom
514 END IF
515
516 END DO ! imol
517
518 END SUBROUTINE get_molecule_set_info
519
520END MODULE molecule_types
Initialize the collective variables types.
recursive subroutine, public colvar_release(colvar)
releases the memory that might have been allocated by the colvar
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Define the molecule kind structure types and the corresponding functionality.
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.
Define the data structure for the molecule information.
subroutine, public deallocate_molecule_set(molecule_set)
Deallocate a molecule set.
subroutine, public deallocate_global_constraint(gci)
Deallocate a global constraint.
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.
subroutine, public allocate_molecule_set(molecule_set, nmolecule)
Allocate a molecule set.
subroutine, public molecule_of_atom(molecule_set, atom_to_mol)
finds for each atom the molecule it belongs to
subroutine, public set_molecule_set(molecule_set, first_atom, last_atom)
Set a molecule data set.
subroutine, public set_molecule(molecule, molecule_kind, lmi, lci, lcolv, lg3x3, lg4x6)
Set a molecule data set.
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.
parameters for a collective variable