LCOV - code coverage report
Current view: top level - src - constraint_fxd.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:4dc10b3) Lines: 207 241 85.9 %
Date: 2024-11-21 06:45:46 Functions: 6 6 100.0 %

          Line data    Source code
       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             : !> \par History
      10             : !>      none
      11             : ! **************************************************************************************************
      12             : MODULE constraint_fxd
      13             : 
      14             :    USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
      15             :    USE atomic_kind_types,               ONLY: get_atomic_kind_set
      16             :    USE cell_types,                      ONLY: use_perd_x,&
      17             :                                               use_perd_xy,&
      18             :                                               use_perd_xyz,&
      19             :                                               use_perd_xz,&
      20             :                                               use_perd_y,&
      21             :                                               use_perd_yz,&
      22             :                                               use_perd_z
      23             :    USE colvar_types,                    ONLY: colvar_type
      24             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      25             :                                               cp_subsys_type
      26             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      27             :    USE force_env_types,                 ONLY: force_env_get,&
      28             :                                               force_env_type
      29             :    USE kinds,                           ONLY: dp
      30             :    USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
      31             :    USE molecule_kind_types,             ONLY: fixd_constraint_type,&
      32             :                                               get_molecule_kind,&
      33             :                                               local_fixd_constraint_type,&
      34             :                                               molecule_kind_type
      35             :    USE molecule_types,                  ONLY: local_g3x3_constraint_type,&
      36             :                                               local_g4x6_constraint_type
      37             :    USE particle_list_types,             ONLY: particle_list_type
      38             :    USE particle_types,                  ONLY: particle_type
      39             :    USE util,                            ONLY: sort
      40             : #include "./base/base_uses.f90"
      41             : 
      42             :    IMPLICIT NONE
      43             : 
      44             :    PRIVATE
      45             :    PUBLIC :: fix_atom_control, &
      46             :              check_fixed_atom_cns_g3x3, &
      47             :              check_fixed_atom_cns_g4x6, &
      48             :              check_fixed_atom_cns_colv, &
      49             :              create_local_fixd_list, &
      50             :              release_local_fixd_list
      51             : 
      52             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_fxd'
      53             : 
      54             : CONTAINS
      55             : 
      56             : ! **************************************************************************************************
      57             : !> \brief allows for fix atom constraints
      58             : !> \param force_env ...
      59             : !> \param w ...
      60             : !> \par History
      61             : !>      - optionally apply fix atom constraint to random forces (Langevin)
      62             : !>        (04.10.206,MK)
      63             : ! **************************************************************************************************
      64      113146 :    SUBROUTINE fix_atom_control(force_env, w)
      65             :       TYPE(force_env_type), POINTER                      :: force_env
      66             :       REAL(KIND=dp), DIMENSION(:, :), OPTIONAL           :: w
      67             : 
      68             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'fix_atom_control'
      69             : 
      70             :       INTEGER :: handle, i, ifixd, ii, ikind, iparticle, iparticle_local, my_atm_fixed, natom, &
      71             :          ncore, nfixed_atoms, nkind, nparticle, nparticle_local, nshell, shell_index
      72             :       LOGICAL                                            :: shell_present
      73      113146 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: force
      74             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
      75             :       TYPE(cp_subsys_type), POINTER                      :: subsys
      76             :       TYPE(distribution_1d_type), POINTER                :: local_particles
      77      113146 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
      78      113146 :       TYPE(local_fixd_constraint_type), POINTER          :: lfixd_list(:)
      79             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
      80      113146 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
      81             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
      82             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
      83             :                                                             shell_particles
      84      113146 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
      85      113146 :                                                             shell_particle_set
      86             : 
      87      113146 :       CALL timeset(routineN, handle)
      88             : 
      89      113146 :       NULLIFY (atomic_kinds)
      90      113146 :       NULLIFY (core_particles)
      91      113146 :       NULLIFY (particles)
      92      113146 :       NULLIFY (shell_particles)
      93             :       shell_present = .FALSE.
      94             : 
      95      113146 :       NULLIFY (lfixd_list)
      96             :       CALL force_env_get(force_env=force_env, &
      97      113146 :                          subsys=subsys)
      98             :       CALL cp_subsys_get(subsys=subsys, &
      99             :                          atomic_kinds=atomic_kinds, &
     100             :                          core_particles=core_particles, &
     101             :                          local_particles=local_particles, &
     102             :                          molecule_kinds=molecule_kinds, &
     103             :                          natom=natom, &
     104             :                          ncore=ncore, &
     105             :                          nshell=nshell, &
     106             :                          particles=particles, &
     107      113146 :                          shell_particles=shell_particles)
     108             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, &
     109      113146 :                                shell_present=shell_present)
     110             : 
     111      113146 :       particle_set => particles%els
     112      113146 :       CPASSERT((SIZE(particle_set) == natom))
     113      113146 :       IF (shell_present) THEN
     114        9446 :          core_particle_set => core_particles%els
     115        9446 :          CPASSERT((SIZE(core_particle_set) == ncore))
     116        9446 :          shell_particle_set => shell_particles%els
     117        9446 :          CPASSERT((SIZE(shell_particle_set) == nshell))
     118             :       END IF
     119      113146 :       nparticle = natom + nshell
     120      113146 :       molecule_kind_set => molecule_kinds%els
     121             : 
     122      113146 :       nkind = molecule_kinds%n_els
     123      113146 :       my_atm_fixed = 0
     124     1963375 :       DO ikind = 1, nkind
     125     1850229 :          molecule_kind => molecule_kind_set(ikind)
     126     1850229 :          CALL get_molecule_kind(molecule_kind, nfixd=nfixed_atoms)
     127     1963375 :          my_atm_fixed = my_atm_fixed + nfixed_atoms
     128             :       END DO
     129             : 
     130      113146 :       IF (my_atm_fixed /= 0) THEN
     131        1226 :          IF (.NOT. PRESENT(w)) THEN
     132             :             ! Allocate scratch array
     133        3630 :             ALLOCATE (force(3, nparticle))
     134     1050354 :             force(:, :) = 0.0_dp
     135        4584 :             DO i = 1, SIZE(local_particles%n_el)
     136        3374 :                nparticle_local = local_particles%n_el(i)
     137      110419 :                DO iparticle_local = 1, nparticle_local
     138      105835 :                   iparticle = local_particles%list(i)%array(iparticle_local)
     139      105835 :                   shell_index = particle_set(iparticle)%shell_index
     140      109209 :                   IF (shell_index == 0) THEN
     141      321964 :                      force(:, iparticle) = particle_set(iparticle)%f(:)
     142             :                   ELSE
     143      101376 :                      force(:, iparticle) = core_particle_set(shell_index)%f(:)
     144      101376 :                      force(:, natom + shell_index) = shell_particle_set(shell_index)%f(:)
     145             :                   END IF
     146             :                END DO
     147             :             END DO
     148             :          END IF
     149             : 
     150             :          ! Create the list of locally fixed atoms
     151        1226 :          CALL create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, local_particles)
     152             : 
     153             :          ! Apply fixed atom constraint
     154       58064 :          DO ifixd = 1, SIZE(lfixd_list)
     155       56838 :             ikind = lfixd_list(ifixd)%ikind
     156       56838 :             ii = lfixd_list(ifixd)%ifixd_index
     157       56838 :             molecule_kind => molecule_kind_set(ikind)
     158       56838 :             CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list)
     159       58064 :             IF (.NOT. fixd_list(ii)%restraint%active) THEN
     160       56508 :                iparticle = fixd_list(ii)%fixd
     161       56508 :                shell_index = particle_set(iparticle)%shell_index
     162             :                ! Select constraint type
     163       56508 :                IF (PRESENT(w)) THEN
     164          39 :                   SELECT CASE (fixd_list(ii)%itype)
     165             :                   CASE (use_perd_x)
     166           0 :                      w(1, iparticle) = 0.0_dp
     167             :                   CASE (use_perd_y)
     168           0 :                      w(2, iparticle) = 0.0_dp
     169             :                   CASE (use_perd_z)
     170           0 :                      w(3, iparticle) = 0.0_dp
     171             :                   CASE (use_perd_xy)
     172           0 :                      w(1, iparticle) = 0.0_dp
     173           0 :                      w(2, iparticle) = 0.0_dp
     174             :                   CASE (use_perd_xz)
     175           0 :                      w(1, iparticle) = 0.0_dp
     176           0 :                      w(3, iparticle) = 0.0_dp
     177             :                   CASE (use_perd_yz)
     178           0 :                      w(2, iparticle) = 0.0_dp
     179           0 :                      w(3, iparticle) = 0.0_dp
     180             :                   CASE (use_perd_xyz)
     181         156 :                      w(:, iparticle) = 0.0_dp
     182             :                   END SELECT
     183             :                ELSE
     184       61845 :                   SELECT CASE (fixd_list(ii)%itype)
     185             :                   CASE (use_perd_x)
     186        5376 :                      force(1, iparticle) = 0.0_dp
     187        5376 :                      IF (shell_index /= 0) THEN
     188           0 :                         force(1, natom + shell_index) = 0.0_dp
     189             :                      END IF
     190             :                   CASE (use_perd_y)
     191        5376 :                      force(2, iparticle) = 0.0_dp
     192        5376 :                      IF (shell_index /= 0) THEN
     193           0 :                         force(2, natom + shell_index) = 0.0_dp
     194             :                      END IF
     195             :                   CASE (use_perd_z)
     196        5376 :                      force(3, iparticle) = 0.0_dp
     197        5376 :                      IF (shell_index /= 0) THEN
     198           0 :                         force(3, natom + shell_index) = 0.0_dp
     199             :                      END IF
     200             :                   CASE (use_perd_xy)
     201        5376 :                      force(1, iparticle) = 0.0_dp
     202        5376 :                      force(2, iparticle) = 0.0_dp
     203        5376 :                      IF (shell_index /= 0) THEN
     204           0 :                         force(1, natom + shell_index) = 0.0_dp
     205           0 :                         force(2, natom + shell_index) = 0.0_dp
     206             :                      END IF
     207             :                   CASE (use_perd_xz)
     208           0 :                      force(1, iparticle) = 0.0_dp
     209           0 :                      force(3, iparticle) = 0.0_dp
     210           0 :                      IF (shell_index /= 0) THEN
     211           0 :                         force(1, natom + shell_index) = 0.0_dp
     212           0 :                         force(3, natom + shell_index) = 0.0_dp
     213             :                      END IF
     214             :                   CASE (use_perd_yz)
     215           0 :                      force(2, iparticle) = 0.0_dp
     216           0 :                      force(3, iparticle) = 0.0_dp
     217           0 :                      IF (shell_index /= 0) THEN
     218           0 :                         force(2, natom + shell_index) = 0.0_dp
     219           0 :                         force(3, natom + shell_index) = 0.0_dp
     220             :                      END IF
     221             :                   CASE (use_perd_xyz)
     222      139860 :                      force(:, iparticle) = 0.0_dp
     223       91434 :                      IF (shell_index /= 0) THEN
     224       88704 :                         force(:, natom + shell_index) = 0.0_dp
     225             :                      END IF
     226             :                   END SELECT
     227             :                END IF
     228             :             END IF
     229             :          END DO
     230        1226 :          CALL release_local_fixd_list(lfixd_list)
     231             : 
     232        1226 :          IF (.NOT. PRESENT(w)) THEN
     233        1210 :             CALL force_env%para_env%sum(force)
     234      212808 :             DO iparticle = 1, natom
     235      211598 :                shell_index = particle_set(iparticle)%shell_index
     236      212808 :                IF (shell_index == 0) THEN
     237      643640 :                   particle_set(iparticle)%f(:) = force(:, iparticle)
     238             :                ELSE
     239      202752 :                   core_particle_set(shell_index)%f(:) = force(:, iparticle)
     240      202752 :                   shell_particle_set(shell_index)%f(:) = force(:, natom + shell_index)
     241             :                END IF
     242             :             END DO
     243        1210 :             DEALLOCATE (force)
     244             :          END IF
     245             :       END IF
     246             : 
     247      113146 :       CALL timestop(handle)
     248             : 
     249      226292 :    END SUBROUTINE fix_atom_control
     250             : 
     251             : ! **************************************************************************************************
     252             : !> \brief ...
     253             : !> \param imass1 ...
     254             : !> \param imass2 ...
     255             : !> \param imass3 ...
     256             : !> \param index_a ...
     257             : !> \param index_b ...
     258             : !> \param index_c ...
     259             : !> \param fixd_list ...
     260             : !> \param lg3x3 ...
     261             : !> \par History
     262             : !>      none
     263             : ! **************************************************************************************************
     264      498309 :    SUBROUTINE check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
     265             :                                         index_a, index_b, index_c, fixd_list, lg3x3)
     266             :       REAL(KIND=dp), INTENT(INOUT)                       :: imass1, imass2, imass3
     267             :       INTEGER, INTENT(IN)                                :: index_a, index_b, index_c
     268             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     269             :       TYPE(local_g3x3_constraint_type)                   :: lg3x3
     270             : 
     271             :       INTEGER                                            :: i
     272             : 
     273      498309 :       IF (lg3x3%init) THEN
     274      485895 :          imass1 = lg3x3%imass1
     275      485895 :          imass2 = lg3x3%imass2
     276      485895 :          imass3 = lg3x3%imass3
     277             :       ELSE
     278       12414 :          IF (ASSOCIATED(fixd_list)) THEN
     279        2317 :             IF (SIZE(fixd_list) > 0) THEN
     280           3 :                DO i = 1, SIZE(fixd_list)
     281           3 :                   IF (fixd_list(i)%fixd == index_a) THEN
     282           2 :                      IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE
     283           2 :                      IF (.NOT. fixd_list(i)%restraint%active) imass1 = 0.0_dp
     284             :                      EXIT
     285             :                   END IF
     286             :                END DO
     287           6 :                DO i = 1, SIZE(fixd_list)
     288           6 :                   IF (fixd_list(i)%fixd == index_b) THEN
     289           0 :                      IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE
     290           0 :                      IF (.NOT. fixd_list(i)%restraint%active) imass2 = 0.0_dp
     291             :                      EXIT
     292             :                   END IF
     293             :                END DO
     294           6 :                DO i = 1, SIZE(fixd_list)
     295           6 :                   IF (fixd_list(i)%fixd == index_c) THEN
     296           0 :                      IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE
     297           0 :                      IF (.NOT. fixd_list(i)%restraint%active) imass3 = 0.0_dp
     298             :                      EXIT
     299             :                   END IF
     300             :                END DO
     301             :             END IF
     302             :          END IF
     303       12414 :          lg3x3%imass1 = imass1
     304       12414 :          lg3x3%imass2 = imass2
     305       12414 :          lg3x3%imass3 = imass3
     306       12414 :          lg3x3%init = .TRUE.
     307             :       END IF
     308      498309 :    END SUBROUTINE check_fixed_atom_cns_g3x3
     309             : 
     310             : ! **************************************************************************************************
     311             : !> \brief ...
     312             : !> \param imass1 ...
     313             : !> \param imass2 ...
     314             : !> \param imass3 ...
     315             : !> \param imass4 ...
     316             : !> \param index_a ...
     317             : !> \param index_b ...
     318             : !> \param index_c ...
     319             : !> \param index_d ...
     320             : !> \param fixd_list ...
     321             : !> \param lg4x6 ...
     322             : !> \par History
     323             : !>      none
     324             : ! **************************************************************************************************
     325        3496 :    SUBROUTINE check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
     326             :                                         index_a, index_b, index_c, index_d, fixd_list, lg4x6)
     327             :       REAL(KIND=dp), INTENT(INOUT)                       :: imass1, imass2, imass3, imass4
     328             :       INTEGER, INTENT(IN)                                :: index_a, index_b, index_c, index_d
     329             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     330             :       TYPE(local_g4x6_constraint_type)                   :: lg4x6
     331             : 
     332             :       INTEGER                                            :: i
     333             : 
     334        3496 :       IF (lg4x6%init) THEN
     335        3300 :          imass1 = lg4x6%imass1
     336        3300 :          imass2 = lg4x6%imass2
     337        3300 :          imass3 = lg4x6%imass3
     338        3300 :          imass4 = lg4x6%imass4
     339             :       ELSE
     340         196 :          IF (ASSOCIATED(fixd_list)) THEN
     341          64 :             IF (SIZE(fixd_list) > 0) THEN
     342        2080 :                DO i = 1, SIZE(fixd_list)
     343        2080 :                   IF (fixd_list(i)%fixd == index_a) THEN
     344          64 :                      IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE
     345          64 :                      IF (.NOT. fixd_list(i)%restraint%active) imass1 = 0.0_dp
     346             :                      EXIT
     347             :                   END IF
     348             :                END DO
     349        4160 :                DO i = 1, SIZE(fixd_list)
     350        4160 :                   IF (fixd_list(i)%fixd == index_b) THEN
     351           0 :                      IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE
     352           0 :                      IF (.NOT. fixd_list(i)%restraint%active) imass2 = 0.0_dp
     353             :                      EXIT
     354             :                   END IF
     355             :                END DO
     356        4160 :                DO i = 1, SIZE(fixd_list)
     357        4160 :                   IF (fixd_list(i)%fixd == index_c) THEN
     358           0 :                      IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE
     359           0 :                      IF (.NOT. fixd_list(i)%restraint%active) imass3 = 0.0_dp
     360             :                      EXIT
     361             :                   END IF
     362             :                END DO
     363        4160 :                DO i = 1, SIZE(fixd_list)
     364        4160 :                   IF (fixd_list(i)%fixd == index_d) THEN
     365           0 :                      IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE
     366           0 :                      IF (.NOT. fixd_list(i)%restraint%active) imass4 = 0.0_dp
     367             :                      EXIT
     368             :                   END IF
     369             :                END DO
     370             :             END IF
     371             :          END IF
     372         196 :          lg4x6%imass1 = imass1
     373         196 :          lg4x6%imass2 = imass2
     374         196 :          lg4x6%imass3 = imass3
     375         196 :          lg4x6%imass4 = imass4
     376         196 :          lg4x6%init = .TRUE.
     377             :       END IF
     378        3496 :    END SUBROUTINE check_fixed_atom_cns_g4x6
     379             : 
     380             : ! **************************************************************************************************
     381             : !> \brief ...
     382             : !> \param fixd_list ...
     383             : !> \param colvar ...
     384             : !> \par History
     385             : !>      none
     386             : ! **************************************************************************************************
     387      404199 :    SUBROUTINE check_fixed_atom_cns_colv(fixd_list, colvar)
     388             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     389             :       TYPE(colvar_type), POINTER                         :: colvar
     390             : 
     391             :       INTEGER                                            :: i, j, k
     392             : 
     393      404199 :       IF (ASSOCIATED(fixd_list)) THEN
     394             :          IF (ASSOCIATED(fixd_list)) THEN
     395         538 :             IF (SIZE(fixd_list) > 0) THEN
     396        2152 :                DO i = 1, SIZE(colvar%i_atom)
     397        1614 :                   j = colvar%i_atom(i)
     398        4552 :                   DO k = 1, SIZE(fixd_list)
     399        4014 :                      IF (fixd_list(k)%fixd == j) THEN
     400         538 :                         IF (fixd_list(k)%itype /= use_perd_xyz) CYCLE
     401         538 :                         IF (.NOT. fixd_list(k)%restraint%active) &
     402        1984 :                            colvar%dsdr(:, i) = 0.0_dp
     403             :                         EXIT
     404             :                      END IF
     405             :                   END DO
     406             :                END DO
     407             :             END IF
     408             :          END IF
     409             :       END IF
     410             : 
     411      404199 :    END SUBROUTINE check_fixed_atom_cns_colv
     412             : 
     413             : ! **************************************************************************************************
     414             : !> \brief setup a list of local atoms on which to apply constraints/restraints
     415             : !> \param lfixd_list ...
     416             : !> \param nkind ...
     417             : !> \param molecule_kind_set ...
     418             : !> \param local_particles ...
     419             : !> \author Teodoro Laino [tlaino] - 11.2008
     420             : ! **************************************************************************************************
     421       99258 :    SUBROUTINE create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, &
     422             :                                      local_particles)
     423             :       TYPE(local_fixd_constraint_type), POINTER          :: lfixd_list(:)
     424             :       INTEGER, INTENT(IN)                                :: nkind
     425             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     426             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     427             : 
     428             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'create_local_fixd_list'
     429             : 
     430             :       INTEGER                                            :: handle, i, ikind, iparticle, &
     431             :                                                             iparticle_local, isize, jsize, ncnst, &
     432             :                                                             nparticle_local, nparticle_local_all, &
     433             :                                                             nsize
     434       99258 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: fixed_atom_all, kind_index_all, &
     435       99258 :                                                             local_particle_all, work0, work1, work2
     436       99258 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     437             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     438             : 
     439       99258 :       CALL timeset(routineN, handle)
     440       99258 :       CPASSERT(.NOT. ASSOCIATED(lfixd_list))
     441       99258 :       nsize = 0
     442     1872453 :       DO ikind = 1, nkind
     443     1773195 :          molecule_kind => molecule_kind_set(ikind)
     444     1773195 :          CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list)
     445     1872453 :          IF (ASSOCIATED(fixd_list)) THEN
     446       82118 :             nsize = nsize + SIZE(fixd_list)
     447             :          END IF
     448             :       END DO
     449       99258 :       IF (nsize /= 0) THEN
     450        7368 :          ALLOCATE (fixed_atom_all(nsize))
     451        4912 :          ALLOCATE (work0(nsize))
     452        4912 :          ALLOCATE (work1(nsize))
     453        4912 :          ALLOCATE (kind_index_all(nsize))
     454        2456 :          nsize = 0
     455       84574 :          DO ikind = 1, nkind
     456       82118 :             molecule_kind => molecule_kind_set(ikind)
     457       82118 :             CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list)
     458       84574 :             IF (ASSOCIATED(fixd_list)) THEN
     459      309708 :                DO i = 1, SIZE(fixd_list)
     460      227590 :                   nsize = nsize + 1
     461      227590 :                   work0(nsize) = i
     462      227590 :                   kind_index_all(nsize) = ikind
     463      309708 :                   fixed_atom_all(nsize) = fixd_list(i)%fixd
     464             :                END DO
     465             :             END IF
     466             :          END DO
     467             :          ! Sort the number of all atoms to be constrained/restrained
     468        2456 :          CALL sort(fixed_atom_all, nsize, work1)
     469             : 
     470             :          ! Sort the local particles
     471        2456 :          nparticle_local_all = 0
     472        9320 :          DO i = 1, SIZE(local_particles%n_el)
     473        9320 :             nparticle_local_all = nparticle_local_all + local_particles%n_el(i)
     474             :          END DO
     475        7190 :          ALLOCATE (local_particle_all(nparticle_local_all))
     476        4734 :          ALLOCATE (work2(nparticle_local_all))
     477        2456 :          nparticle_local_all = 0
     478        9320 :          DO i = 1, SIZE(local_particles%n_el)
     479        6864 :             nparticle_local = local_particles%n_el(i)
     480      244037 :             DO iparticle_local = 1, nparticle_local
     481      234717 :                nparticle_local_all = nparticle_local_all + 1
     482      234717 :                iparticle = local_particles%list(i)%array(iparticle_local)
     483      241581 :                local_particle_all(nparticle_local_all) = iparticle
     484             :             END DO
     485             :          END DO
     486        2456 :          CALL sort(local_particle_all, nparticle_local_all, work2)
     487             : 
     488             :          ! Count the amount of local constraints/restraints
     489        2456 :          ncnst = 0
     490        2456 :          jsize = 1
     491      174275 :          Loop_count: DO isize = 1, nparticle_local_all
     492      396079 :             DO WHILE (local_particle_all(isize) > fixed_atom_all(jsize))
     493      224260 :                jsize = jsize + 1
     494      396079 :                IF (jsize > nsize) THEN
     495        2456 :                   jsize = nsize
     496             :                   EXIT Loop_count
     497             :                END IF
     498             :             END DO
     499      173186 :             IF (local_particle_all(isize) == fixed_atom_all(jsize)) ncnst = ncnst + 1
     500             :          END DO Loop_count
     501             : 
     502             :          ! Allocate local fixed atom array
     503      120810 :          ALLOCATE (lfixd_list(ncnst))
     504             : 
     505             :          ! Fill array with constraints infos
     506             :          ncnst = 0
     507             :          jsize = 1
     508      174275 :          Loop_fill: DO isize = 1, nparticle_local_all
     509      396079 :             DO WHILE (local_particle_all(isize) > fixed_atom_all(jsize))
     510      224260 :                jsize = jsize + 1
     511      396079 :                IF (jsize > nsize) THEN
     512        2456 :                   jsize = nsize
     513             :                   EXIT Loop_fill
     514             :                END IF
     515             :             END DO
     516      173186 :             IF (local_particle_all(isize) == fixed_atom_all(jsize)) THEN
     517      113829 :                ncnst = ncnst + 1
     518      113829 :                lfixd_list(ncnst)%ifixd_index = work0(work1(jsize))
     519      113829 :                lfixd_list(ncnst)%ikind = kind_index_all(work1(jsize))
     520             :             END IF
     521             :          END DO Loop_fill
     522             : 
     523             :          ! Deallocate working arrays
     524        2456 :          DEALLOCATE (local_particle_all)
     525        2456 :          DEALLOCATE (work2)
     526        2456 :          DEALLOCATE (fixed_atom_all)
     527        2456 :          DEALLOCATE (work1)
     528        2456 :          DEALLOCATE (kind_index_all)
     529             :       ELSE
     530             :          ! Allocate local fixed atom array with dimension 0
     531       96802 :          ALLOCATE (lfixd_list(0))
     532             :       END IF
     533       99258 :       CALL timestop(handle)
     534      198516 :    END SUBROUTINE create_local_fixd_list
     535             : 
     536             : ! **************************************************************************************************
     537             : !> \brief destroy the list of local atoms on which to apply constraints/restraints
     538             : !>      Teodoro Laino [tlaino] - 11.2008
     539             : !> \param lfixd_list ...
     540             : ! **************************************************************************************************
     541       99258 :    SUBROUTINE release_local_fixd_list(lfixd_list)
     542             :       TYPE(local_fixd_constraint_type), POINTER          :: lfixd_list(:)
     543             : 
     544       99258 :       CPASSERT(ASSOCIATED(lfixd_list))
     545       99258 :       DEALLOCATE (lfixd_list)
     546       99258 :    END SUBROUTINE release_local_fixd_list
     547             : 
     548             : END MODULE constraint_fxd

Generated by: LCOV version 1.15