CP2K 2.4 (Revision 12889)

fist_neighbor_list_control.f90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------------!
00002 !   CP2K: A general program to perform molecular dynamics simulations         !
00003 !   Copyright (C) 2000 - 2013  CP2K developers group                          !
00004 !-----------------------------------------------------------------------------!
00005 
00006 ! *****************************************************************************
00013 MODULE fist_neighbor_list_control
00014 
00015   USE atomic_kind_types,               ONLY: atomic_kind_type,&
00016                                              get_atomic_kind_set
00017   USE cell_types,                      ONLY: cell_clone,&
00018                                              cell_create,&
00019                                              cell_release,&
00020                                              cell_type
00021   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
00022                                              cp_print_key_unit_nr
00023   USE cp_para_types,                   ONLY: cp_para_env_type
00024   USE distribution_1d_types,           ONLY: distribution_1d_type
00025   USE exclusion_types,                 ONLY: exclusion_type
00026   USE f77_blas
00027   USE fist_neighbor_list_types,        ONLY: fist_neighbor_type
00028   USE fist_neighbor_lists,             ONLY: build_fist_neighbor_lists
00029   USE fist_nonbond_env_types,          ONLY: fist_nonbond_env_get,&
00030                                              fist_nonbond_env_set,&
00031                                              fist_nonbond_env_type,&
00032                                              pos_type
00033   USE input_section_types,             ONLY: section_vals_type,&
00034                                              section_vals_val_get
00035   USE kinds,                           ONLY: dp
00036   USE message_passing,                 ONLY: mp_max
00037   USE pair_potential_types,            ONLY: pair_potential_pp_type,&
00038                                              siepmann_type,&
00039                                              tersoff_type
00040   USE particle_types,                  ONLY: particle_type
00041   USE timings,                         ONLY: timeset,&
00042                                              timestop
00043 #include "cp_common_uses.h"
00044 
00045   IMPLICIT NONE
00046 
00047   PRIVATE
00048 
00049   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_neighbor_list_control'
00050 
00051   PUBLIC :: list_control
00052 
00053 !***
00054 
00055 CONTAINS
00056 
00057 ! to decide whether the neighbor list is to be updated or not
00058 ! based on a displacement criterion;
00059 ! if any particle has moved by 0.5*verlet_skin from the previous
00060 ! list update, then the list routine is called.
00061 
00062 ! *****************************************************************************
00063   SUBROUTINE list_control ( atomic_kind_set, particle_set, local_particles, &
00064        cell, fist_nonbond_env, para_env, mm_section,shell_particle_set,&
00065        core_particle_set, force_update, exclusions, error)
00066 
00067     TYPE(atomic_kind_type), POINTER          :: atomic_kind_set(:)
00068     TYPE(particle_type), POINTER             :: particle_set(:)
00069     TYPE(distribution_1d_type), POINTER      :: local_particles
00070     TYPE(cell_type), POINTER                 :: cell
00071     TYPE(fist_nonbond_env_type), POINTER     :: fist_nonbond_env
00072     TYPE(cp_para_env_type), POINTER          :: para_env
00073     TYPE(section_vals_type), POINTER         :: mm_section
00074     TYPE(particle_type), OPTIONAL, POINTER   :: shell_particle_set(:), 
00075                                                 core_particle_set(:)
00076     LOGICAL, INTENT(IN), OPTIONAL            :: force_update
00077     TYPE(exclusion_type), DIMENSION(:), 
00078       OPTIONAL                               :: exclusions
00079     TYPE(cp_error_type), INTENT(inout)       :: error
00080 
00081     CHARACTER(LEN=*), PARAMETER :: routineN = 'list_control', 
00082       routineP = moduleN//':'//routineN
00083 
00084     INTEGER :: counter, handle, ikind, iparticle, iparticle_kind, 
00085       iparticle_local, ishell, jkind, last_update, nparticle, nparticle_kind, 
00086       nparticle_local, nshell, num_update, output_unit, stat
00087     LOGICAL :: build_from_scratch, failure, geo_check, shell_adiabatic, 
00088       shell_present, update_neighbor_lists
00089     LOGICAL, DIMENSION(:, :), POINTER        :: full_nl
00090     REAL(KIND=dp)                            :: aup, dr2, dr2_max, 
00091                                                 ei_scale14, lup, vdw_scale14, 
00092                                                 verlet_skin
00093     REAL(KIND=dp), DIMENSION(3)              :: dr, rab, s, s2r
00094     REAL(KIND=dp), DIMENSION(:, :), POINTER  :: rlist_cut, rlist_lowsq
00095     TYPE(cell_type), POINTER                 :: cell_last_update
00096     TYPE(cp_logger_type), POINTER            :: logger
00097     TYPE(fist_neighbor_type), POINTER        :: nonbonded
00098     TYPE(pair_potential_pp_type), POINTER    :: potparm
00099     TYPE(pos_type), DIMENSION(:), POINTER    :: r_last_update, 
00100                                                 r_last_update_pbc, 
00101                                                 rcore_last_update_pbc, 
00102                                                 rshell_last_update_pbc
00103 
00104     failure = .FALSE.
00105     CALL timeset(routineN,handle)
00106     NULLIFY(logger)
00107     logger => cp_error_get_logger(error)
00108 
00109     ! *** Assigning local pointers ***
00110     CALL fist_nonbond_env_get(fist_nonbond_env,&
00111                               nonbonded=nonbonded,&
00112                               rlist_cut=rlist_cut,&
00113                               rlist_lowsq=rlist_lowsq,&
00114                               aup=aup,&
00115                               lup=lup, &
00116                               ei_scale14=ei_scale14, &
00117                               vdw_scale14=vdw_scale14, &
00118                               counter=counter,&
00119                               r_last_update=r_last_update,&
00120                               r_last_update_pbc=r_last_update_pbc,&
00121                               rshell_last_update_pbc=rshell_last_update_pbc,&
00122                               rcore_last_update_pbc=rcore_last_update_pbc,&
00123                               cell_last_update=cell_last_update,&
00124                               num_update=num_update,&
00125                               potparm=potparm,&
00126                               last_update=last_update,error=error)
00127 
00128     nparticle      = SIZE(particle_set)
00129     nparticle_kind = SIZE(atomic_kind_set)
00130     nshell = 0
00131     CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,&
00132          shell_present=shell_present, shell_adiabatic=shell_adiabatic)
00133     IF(shell_present) THEN
00134        nshell = SIZE(shell_particle_set)
00135     END IF
00136 
00137     ! *** Check, if the neighbor lists have to be built or updated ***
00138     update_neighbor_lists = .FALSE.
00139     CALL section_vals_val_get(mm_section,"NEIGHBOR_LISTS%NEIGHBOR_LISTS_FROM_SCRATCH",&
00140          l_val=build_from_scratch,error=error)
00141     CALL section_vals_val_get(mm_section,"NEIGHBOR_LISTS%GEO_CHECK",&
00142          l_val=geo_check,error=error)
00143     IF (ASSOCIATED(r_last_update)) THEN
00144        ! Determine the maximum of the squared displacement, compared to
00145        ! r_last_update.
00146        CALL section_vals_val_get(mm_section,"NEIGHBOR_LISTS%VERLET_SKIN",&
00147             r_val=verlet_skin,error=error)
00148        dr2_max = 0.0_dp
00149        DO iparticle_kind=1,nparticle_kind
00150           nparticle_local = local_particles%n_el(iparticle_kind)
00151           DO iparticle_local=1,nparticle_local
00152              iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
00153              s2r = r_last_update(iparticle)%r
00154              s=particle_set(iparticle)%r(:)
00155              dr(:) = s2r - s
00156              dr2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
00157              dr2_max = MAX(dr2_max,dr2)
00158           END DO
00159        END DO
00160 
00161        CALL mp_max(dr2_max,para_env%group)
00162 
00163        ! If the maximum distplacement is too large, ...
00164        IF (dr2_max > 0.25_dp*verlet_skin**2.OR.build_from_scratch) THEN
00165           DO iparticle=1,nparticle
00166              r_last_update(iparticle )%r=particle_set(iparticle)%r(:)
00167           END DO
00168           update_neighbor_lists = .TRUE.
00169        END IF
00170     ELSE
00171        ! There is no r_last_update to compare with. Neighbor lists from scratch.
00172        ALLOCATE (r_last_update(nparticle),STAT=stat)
00173        CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00174        DO iparticle=1,nparticle
00175           r_last_update(iparticle )%r=particle_set(iparticle)%r(:)
00176        END DO
00177 
00178        update_neighbor_lists = .TRUE.
00179        build_from_scratch = .TRUE.
00180     END IF
00181     ! Force Update
00182     IF (PRESENT(force_update)) THEN
00183        IF (force_update) update_neighbor_lists = .TRUE.
00184     END IF
00185 
00186     ! Allocate the r_last_update_pbc, rshell_last_update_pbc, rcore_last_update_pbc
00187     IF (.NOT.ASSOCIATED(r_last_update_pbc)) THEN
00188        ALLOCATE (r_last_update_pbc(nparticle),STAT=stat)
00189        CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00190     END IF
00191     IF(shell_present .AND. .NOT.ASSOCIATED(rshell_last_update_pbc)) THEN
00192        ALLOCATE (rshell_last_update_pbc(nshell),STAT=stat)
00193        CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00194     END IF
00195     IF(shell_present .AND. .NOT.ASSOCIATED(rcore_last_update_pbc)) THEN
00196        ALLOCATE (rcore_last_update_pbc(nshell),STAT=stat)
00197        CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00198     END IF
00199 
00200     ! update the neighbor lists
00201     IF (update_neighbor_lists) THEN
00202        ! determine which pairs of atom kinds need full neighbor lists. Full
00203        ! means that atom a is in the neighbor list of atom b and vice versa.
00204        ALLOCATE(full_nl(nparticle_kind,nparticle_kind),stat=stat)
00205        CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00206        IF (ASSOCIATED(potparm)) THEN
00207           DO ikind = 1, nparticle_kind
00208              DO jkind = ikind, nparticle_kind
00209                 full_nl(ikind, jkind) = .FALSE.
00210                 IF (ANY(potparm%pot(ikind, jkind)%pot%type==tersoff_type)) THEN
00211                    full_nl(ikind, jkind) = .TRUE.
00212                 END IF
00213                 IF (ANY(potparm%pot(ikind, jkind)%pot%type==siepmann_type)) THEN
00214                    full_nl(ikind, jkind) = .TRUE.
00215                 END IF
00216                 full_nl(jkind, ikind) = full_nl(ikind, jkind)
00217              END DO
00218           END DO
00219        ELSE
00220           full_nl = .FALSE.
00221        END IF
00222        CALL build_fist_neighbor_lists(atomic_kind_set, particle_set, &
00223             local_particles, cell, rlist_cut, rlist_lowsq, ei_scale14, &
00224             vdw_scale14, nonbonded, para_env, &
00225             build_from_scratch=build_from_scratch, geo_check=geo_check, &
00226             mm_section=mm_section, full_nl=full_nl,&
00227             exclusions=exclusions, error=error)
00228 
00229        CALL cell_release(cell_last_update,error=error)
00230        CALL cell_create(cell_last_update,error=error)
00231        CALL cell_clone(cell,cell_last_update,error)
00232 
00233        IF ( counter > 0 ) THEN
00234           num_update = num_update + 1
00235           lup = counter + 1 - last_update
00236           last_update = counter + 1
00237           aup = aup + (lup - aup)/REAL(num_update,KIND=dp)
00238        ELSE
00239           num_update = 0
00240           lup = 0
00241           last_update = 1
00242           aup = 0.0_dp
00243        END IF
00244 
00245        CALL fist_nonbond_env_set(fist_nonbond_env,&
00246                                  lup=lup,&
00247                                  aup=aup,&
00248                                  r_last_update=r_last_update,&
00249                                  r_last_update_pbc=r_last_update_pbc,&
00250                                  rshell_last_update_pbc=rshell_last_update_pbc,&
00251                                  rcore_last_update_pbc=rcore_last_update_pbc,&
00252                                  nonbonded=nonbonded,&
00253                                  num_update=num_update,&
00254                                  last_update=last_update,&
00255                                  cell_last_update=cell_last_update,error=error)
00256 
00257        output_unit = cp_print_key_unit_nr(logger,mm_section,"PRINT%NEIGHBOR_LISTS",&
00258             extension=".mmLog",error=error)
00259        IF (output_unit>0) THEN
00260           WRITE (UNIT=output_unit,&
00261                FMT="(/,T2,A,/,T52,A,/,A,T31,A,T49,2(1X,F15.2),/,T2,A,/)")&
00262                REPEAT("*",79),"INSTANTANEOUS        AVERAGES",&
00263                " LIST UPDATES[steps]","= ",lup,aup,REPEAT("*",79)
00264        END IF
00265        CALL cp_print_key_finished_output(output_unit,logger,mm_section,&
00266             "PRINT%NEIGHBOR_LISTS",error=error)
00267        DEALLOCATE(full_nl,stat=stat)
00268        CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00269     END IF
00270 
00271     ! Store particle positions after the last update, translated to the
00272     ! primitive cell, in r_last_update_pbc.
00273     DO iparticle=1,nparticle
00274        rab = r_last_update(iparticle)%r
00275        IF (cell%orthorhombic) THEN
00276           rab(1) = - cell%hmat(1,1)*cell%perd(1)*ANINT(cell_last_update%h_inv(1,1)*rab(1))
00277           rab(2) = - cell%hmat(2,2)*cell%perd(2)*ANINT(cell_last_update%h_inv(2,2)*rab(2))
00278           rab(3) = - cell%hmat(3,3)*cell%perd(3)*ANINT(cell_last_update%h_inv(3,3)*rab(3))
00279        ELSE
00280           s(1) = cell_last_update%h_inv(1,1)*rab(1) + cell_last_update%h_inv(1,2)*rab(2) +&
00281                cell_last_update%h_inv(1,3)*rab(3)
00282           s(2) = cell_last_update%h_inv(2,1)*rab(1) + cell_last_update%h_inv(2,2)*rab(2) +&
00283                cell_last_update%h_inv(2,3)*rab(3)
00284           s(3) = cell_last_update%h_inv(3,1)*rab(1) + cell_last_update%h_inv(3,2)*rab(2) +&
00285                cell_last_update%h_inv(3,3)*rab(3)
00286           s(1) = - cell%perd(1)*ANINT(s(1))
00287           s(2) = - cell%perd(2)*ANINT(s(2))
00288           s(3) = - cell%perd(3)*ANINT(s(3))
00289           rab(1) = + cell%hmat(1,1)*s(1) + cell%hmat(1,2)*s(2) + cell%hmat(1,3)*s(3)
00290           rab(2) = + cell%hmat(2,1)*s(1) + cell%hmat(2,2)*s(2) + cell%hmat(2,3)*s(3)
00291           rab(3) = + cell%hmat(3,1)*s(1) + cell%hmat(3,2)*s(2) + cell%hmat(3,3)*s(3)
00292        END IF
00293        r_last_update_pbc(iparticle )%r=particle_set(iparticle)%r+rab
00294        ! Use the same translation for core and shell.
00295        ishell = particle_set(iparticle)%shell_index
00296        IF(ishell/=0) THEN
00297           rshell_last_update_pbc(ishell )%r =  rab + shell_particle_set(ishell)%r(:)
00298           IF(shell_adiabatic) THEN
00299              rcore_last_update_pbc(ishell )%r = rab + core_particle_set(ishell)%r(:)
00300           ELSE
00301              rcore_last_update_pbc(ishell )%r = r_last_update_pbc(iparticle )%r(:)
00302           END IF
00303        END IF
00304     END DO
00305 
00306     counter = counter + 1
00307     CALL fist_nonbond_env_set(fist_nonbond_env,counter=counter,error=error)
00308     CALL timestop(handle)
00309 
00310   END SUBROUTINE list_control
00311 
00312 END MODULE fist_neighbor_list_control