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