|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00009 MODULE thermostat_mapping 00010 00011 USE cp_para_types, ONLY: cp_para_env_type 00012 USE distribution_1d_types, ONLY: distribution_1d_type 00013 USE extended_system_types, ONLY: map_info_type 00014 USE f77_blas 00015 USE input_constants, ONLY: & 00016 do_region_defined, do_region_global, do_region_massive, & 00017 do_region_molecule, do_thermo_communication, & 00018 do_thermo_no_communication, use_perd_x, use_perd_xy, use_perd_xyz, & 00019 use_perd_xz, use_perd_y, use_perd_yz, use_perd_z 00020 USE kinds, ONLY: dp,& 00021 int_size 00022 USE memory_utilities, ONLY: reallocate 00023 USE message_passing, ONLY: mp_allgather,& 00024 mp_bcast,& 00025 mp_sum 00026 USE molecule_kind_types, ONLY: colvar_constraint_type,& 00027 fixd_constraint_type,& 00028 g3x3_constraint_type,& 00029 g4x6_constraint_type,& 00030 get_molecule_kind,& 00031 molecule_kind_type 00032 USE molecule_types_new, ONLY: get_molecule,& 00033 global_constraint_type,& 00034 molecule_type 00035 USE simpar_types, ONLY: simpar_type 00036 USE termination, ONLY: stop_memory,& 00037 stop_program 00038 USE timings, ONLY: timeset,& 00039 timestop 00040 USE util, ONLY: locate,& 00041 sort 00042 #include "cp_common_uses.h" 00043 00044 IMPLICIT NONE 00045 00046 PUBLIC :: thermostat_mapping_region,& 00047 adiabatic_mapping_region,& 00048 init_baro_map_info 00049 00050 PRIVATE 00051 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermostat_mapping' 00052 00053 CONTAINS 00054 ! ***************************************************************************** 00058 SUBROUTINE adiabatic_mapping_region ( map_info, deg_of_freedom, massive_atom_list,& 00059 molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, simpar, & 00060 number, region, gci, shell, map_loc_thermo_gen, sum_of_thermostats, error) 00061 00062 TYPE(map_info_type), POINTER :: map_info 00063 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, 00064 massive_atom_list 00065 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 00066 TYPE(distribution_1d_type), POINTER :: local_molecules 00067 TYPE(molecule_type), POINTER :: molecule_set(:) 00068 TYPE(cp_para_env_type), POINTER :: para_env 00069 INTEGER, INTENT(OUT) :: natoms_local 00070 TYPE(simpar_type), POINTER :: simpar 00071 INTEGER, INTENT(INOUT) :: number 00072 INTEGER, INTENT(IN) :: region 00073 TYPE(global_constraint_type), POINTER :: gci 00074 LOGICAL, INTENT(IN) :: shell 00075 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen 00076 INTEGER, INTENT(INOUT) :: sum_of_thermostats 00077 TYPE(cp_error_type), INTENT(inout) :: error 00078 00079 CHARACTER(LEN=*), PARAMETER :: routineN = 'adiabatic_mapping_region', 00080 routineP = moduleN//':'//routineN 00081 00082 INTEGER :: handle, nkind, nmol_local, 00083 nsize, number_of_thermostats, 00084 stat 00085 INTEGER, DIMENSION(:), POINTER :: const_mol, tot_const 00086 INTEGER, DIMENSION(:, :), POINTER :: point 00087 LOGICAL :: failure 00088 00089 CALL timeset(routineN,handle) 00090 00091 failure = .FALSE. 00092 NULLIFY ( const_mol, tot_const, point) 00093 CPPostcondition(.NOT.ASSOCIATED(deg_of_freedom),cp_failure_level,routineP,error,failure) 00094 CPPostcondition(.NOT.ASSOCIATED(massive_atom_list),cp_failure_level,routineP,error,failure) 00095 00096 nkind = SIZE(molecule_kind_set) 00097 CALL adiabatic_region_evaluate(map_info%dis_type, natoms_local, nmol_local,& 00098 const_mol, tot_const, point, local_molecules, molecule_kind_set, molecule_set,& 00099 region, simpar, shell, sum_of_thermostats, para_env, error) 00100 00101 ! Now we can allocate the target array s_kin and p_kin.. 00102 SELECT CASE(region) 00103 CASE(do_region_global, do_region_molecule, do_region_massive) 00104 nsize = number 00105 CASE DEFAULT 00106 ! STOP PROGRAM 00107 END SELECT 00108 ALLOCATE (map_info%s_kin(nsize),STAT=stat) 00109 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00110 ALLOCATE (map_info%v_scale(nsize),STAT=stat) 00111 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00112 ALLOCATE (map_info%p_kin(3,natoms_local),STAT=stat) 00113 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00114 ALLOCATE (map_info%p_scale(3,natoms_local),STAT=stat) 00115 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00116 ! nullify thermostat pointers 00117 ! Allocate index array to 1 00118 ALLOCATE ( map_info%index(1), STAT=stat) 00119 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00120 ALLOCATE ( map_info%map_index(1), STAT=stat) 00121 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00122 ALLOCATE (deg_of_freedom(1),STAT=stat) 00123 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00124 00125 CALL massive_list_generate (molecule_set, molecule_kind_set, & 00126 local_molecules, para_env, massive_atom_list, region, shell,& 00127 error) 00128 00129 CALL adiabatic_mapping_region_low(region, map_info, nkind, point,& 00130 deg_of_freedom, local_molecules, const_mol, massive_atom_list,& 00131 tot_const, molecule_set, number_of_thermostats, shell, gci,& 00132 map_loc_thermo_gen, error) 00133 00134 number = number_of_thermostats 00135 sum_of_thermostats=number 00136 CALL mp_sum ( sum_of_thermostats, para_env%group ) 00137 00138 ! check = (number==number_of_thermostats) 00139 ! CPPrecondition(check,cp_fatal_level,routineP,error,failure) 00140 DEALLOCATE (const_mol,STAT=stat) 00141 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00142 DEALLOCATE (tot_const,STAT=stat) 00143 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00144 DEALLOCATE (point, STAT = stat ) 00145 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00146 00147 CALL timestop(handle) 00148 00149 END SUBROUTINE adiabatic_mapping_region 00150 00151 ! ***************************************************************************** 00155 SUBROUTINE adiabatic_mapping_region_low( region, map_info, nkind, point,& 00156 deg_of_freedom, local_molecules, const_mol, massive_atom_list, tot_const,& 00157 molecule_set, ntherm, shell, gci, map_loc_thermo_gen,error) 00158 00159 INTEGER, INTENT(IN) :: region 00160 TYPE(map_info_type), POINTER :: map_info 00161 INTEGER :: nkind 00162 INTEGER, DIMENSION(:, :), POINTER :: point 00163 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom 00164 TYPE(distribution_1d_type), POINTER :: local_molecules 00165 INTEGER, DIMENSION(:), POINTER :: const_mol, massive_atom_list, 00166 tot_const 00167 TYPE(molecule_type), POINTER :: molecule_set(:) 00168 INTEGER, INTENT(OUT) :: ntherm 00169 LOGICAL, INTENT(IN) :: shell 00170 TYPE(global_constraint_type), POINTER :: gci 00171 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen 00172 TYPE(cp_error_type), INTENT(inout) :: error 00173 00174 CHARACTER(LEN=*), PARAMETER :: routineN = 'adiabatic_mapping_region_low', 00175 routineP = moduleN//':'//routineN 00176 00177 INTEGER :: first_atom, first_shell, glob_therm_num, handle, icount, 00178 ielement, ii, ikind, imol, imol_local, ipart, jj, k, kk, last_atom, 00179 last_shell, nglob_cns, nmol_local, number 00180 LOGICAL :: check, failure, 00181 global_constraints, 00182 have_thermostat 00183 REAL(dp), SAVE, TARGET :: unity 00184 TYPE(molecule_type), POINTER :: molecule 00185 00186 CALL timeset(routineN,handle) 00187 unity = 1.0_dp 00188 failure = .FALSE. 00189 global_constraints = ASSOCIATED(gci) 00190 deg_of_freedom = 0 00191 icount = 0 00192 number = 0 00193 ntherm = 0 00194 nglob_cns = 0 00195 IF (global_constraints) nglob_cns = gci%ntot-gci%nrestraint 00196 IF ( region==do_region_global) THEN 00197 ! Global Region 00198 check = ( map_info%dis_type == do_thermo_communication ) 00199 CPPostcondition(check,cp_failure_level,routineP,error,failure) 00200 DO ikind = 1, nkind 00201 DO jj = point ( 1, ikind ), point ( 2, ikind ) 00202 IF ( map_loc_thermo_gen ( jj ) /= HUGE ( 0 ) ) THEN 00203 DO ii = 1, 3 00204 map_info%p_kin(ii,jj)%point => map_info%s_kin(1) 00205 map_info%p_scale(ii,jj)%point => map_info%v_scale(1) 00206 END DO 00207 ELSE 00208 DO ii = 1, 3 00209 NULLIFY ( map_info%p_kin(ii,jj)%point ) 00210 map_info%p_scale(ii,jj)%point => unity 00211 END DO 00212 ENDIF 00213 END DO 00214 deg_of_freedom(1) = deg_of_freedom(1) + tot_const(ikind) 00215 map_info%index(1) = 1 00216 map_info%map_index(1) = 1 00217 number = 1 00218 END DO 00219 deg_of_freedom(1) = deg_of_freedom(1) + nglob_cns 00220 ELSE IF ( region==do_region_molecule) THEN 00221 ! Molecular Region 00222 IF ( map_info%dis_type == do_thermo_no_communication ) THEN 00223 ! This is the standard case.. 00224 DO ikind = 1, nkind 00225 nmol_local = local_molecules % n_el (ikind) 00226 DO imol_local = 1, nmol_local 00227 imol = local_molecules%list(ikind)%array(imol_local) 00228 number = number + 1 00229 have_thermostat=.TRUE. 00230 ! determine if the local molecule belongs to a thermostat 00231 DO kk = point ( 1, number ), point ( 2, number ) 00232 ! WRITE ( *, * ) 'kk', size(map_loc_thermo_gen), kk, map_loc_thermo_gen ( kk ) 00233 IF ( map_loc_thermo_gen ( kk ) == HUGE ( 0 ) ) THEN 00234 have_thermostat=.FALSE. 00235 EXIT 00236 ENDIF 00237 END DO 00238 ! If molecule has a thermostat, then map 00239 IF ( have_thermostat ) THEN 00240 ! glob_therm_num is the global thermostat number attached to the local molecule 00241 ! We can test to make sure all atoms in the local molecule belong to the same 00242 ! global thermostat as a way to detect errors. 00243 glob_therm_num = map_loc_thermo_gen ( point ( 1, number ) ) 00244 ntherm = ntherm + 1 00245 CALL reallocate ( map_info%index, 1, ntherm ) 00246 CALL reallocate ( map_info%map_index, 1, ntherm ) 00247 CALL reallocate ( deg_of_freedom, 1, ntherm ) 00248 map_info%index(ntherm) = glob_therm_num 00249 map_info%map_index(ntherm) = ntherm 00250 deg_of_freedom ( ntherm ) = const_mol ( number ) 00251 DO kk = point ( 1, number ), point ( 2, number ) 00252 DO jj = 1, 3 00253 map_info%p_kin(jj,kk) %point => map_info%s_kin(ntherm) 00254 map_info%p_scale(jj,kk) %point => map_info%v_scale(ntherm) 00255 END DO 00256 END DO 00257 ! If molecule has no thermostat, then nullify pointers to that atom 00258 ELSE 00259 DO kk = point ( 1, number ), point ( 2, number ) 00260 DO jj = 1, 3 00261 NULLIFY ( map_info%p_kin(jj,kk) %point ) 00262 map_info%p_scale(jj,kk) %point => unity 00263 END DO 00264 END DO 00265 END IF 00266 END DO 00267 END DO 00268 ELSE IF ( map_info%dis_type == do_thermo_communication ) THEN 00269 ! This case is quite rare and happens only when we have one molecular 00270 ! kind and one molecule.. 00271 CPPostcondition(nkind==1,cp_failure_level,routineP,error,failure) 00272 number = number + 1 00273 ntherm = ntherm + 1 00274 map_info%index(ntherm) = ntherm 00275 map_info%map_index(ntherm) = ntherm 00276 deg_of_freedom ( ntherm ) = deg_of_freedom ( ntherm ) + tot_const ( nkind ) 00277 DO kk = point ( 1, nkind ), point ( 2, nkind ) 00278 IF ( map_loc_thermo_gen ( kk ) /= HUGE ( 0 ) ) THEN 00279 DO jj = 1, 3 00280 map_info%p_kin ( jj, kk ) % point => map_info%s_kin(ntherm) 00281 map_info%p_scale ( jj, kk ) % point => map_info%v_scale(ntherm) 00282 END DO 00283 ELSE 00284 END IF 00285 DO jj = 1, 3 00286 NULLIFY ( map_info%p_kin ( jj, kk ) % point ) 00287 map_info%p_scale ( jj, kk ) % point => unity 00288 END DO 00289 END DO 00290 ELSE 00291 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 00292 END IF 00293 IF (nglob_cns/=0) THEN 00294 CALL stop_program(routineN,moduleN,__LINE__,& 00295 "Molecular thermostats with global constraints are impossible!") 00296 END IF 00297 ELSE IF ( region==do_region_massive) THEN 00298 ! Massive Region 00299 check = ( map_info%dis_type == do_thermo_no_communication ) 00300 CPPostcondition(check,cp_failure_level,routineP,error,failure) 00301 DO ikind = 1, nkind 00302 nmol_local = local_molecules % n_el ( ikind ) 00303 DO imol_local = 1, nmol_local 00304 icount = icount + 1 00305 imol = local_molecules % list ( ikind ) % array ( imol_local ) 00306 molecule => molecule_set ( imol ) 00307 CALL get_molecule ( molecule, first_atom=first_atom, last_atom=last_atom,& 00308 first_shell=first_shell, last_shell=last_shell) 00309 IF (shell) THEN 00310 first_atom = first_shell 00311 last_atom = last_shell 00312 ELSE 00313 IF ((tot_const(icount)>0).OR.(nglob_cns/=0)) THEN 00314 CALL stop_program(routineN,moduleN,__LINE__,& 00315 "Massive thermostats with constraints are impossible!") 00316 END IF 00317 END IF 00318 k = 0 00319 00320 have_thermostat=.TRUE. 00321 DO ii = point ( 1, icount ), point ( 2, icount ) 00322 IF ( map_loc_thermo_gen ( ii ) /= 1 ) THEN 00323 have_thermostat = .FALSE. 00324 EXIT 00325 END IF 00326 END DO 00327 00328 IF ( have_thermostat ) THEN 00329 DO ii = point ( 1, icount ), point ( 2, icount ) 00330 ipart = first_atom + k 00331 ielement = locate(massive_atom_list,ipart) 00332 k = k + 1 00333 DO jj = 1, 3 00334 ntherm = ntherm + 1 00335 CALL reallocate ( map_info%index, 1, ntherm ) 00336 CALL reallocate ( map_info%map_index, 1, ntherm ) 00337 map_info%index(ntherm) = (ielement - 1)*3 + jj 00338 map_info%map_index(ntherm) = ntherm 00339 map_info%p_kin(jj,ii) %point => map_info%s_kin(ntherm) 00340 map_info%p_scale(jj,ii) %point => map_info%v_scale(ntherm) 00341 END DO 00342 END DO 00343 ELSE 00344 DO ii = point ( 1, icount ), point ( 2, icount ) 00345 DO jj = 1, 3 00346 NULLIFY ( map_info%p_kin(jj,ii) %point ) 00347 map_info%p_scale(jj,ii) %point => unity 00348 END DO 00349 END DO 00350 END IF 00351 IF ( first_atom + k -1 /= last_atom ) THEN 00352 CALL stop_program(routineN,moduleN,__LINE__,& 00353 "Inconsistent mapping of particles") 00354 END IF 00355 END DO 00356 END DO 00357 ELSE 00358 CALL stop_program(routineN,moduleN,__LINE__,"Invalid region!") 00359 END IF 00360 00361 CALL timestop(handle) 00362 00363 END SUBROUTINE adiabatic_mapping_region_low 00364 00365 00366 ! ***************************************************************************** 00370 SUBROUTINE adiabatic_region_evaluate(dis_type, natoms_local, nmol_local, const_mol,& 00371 tot_const, point, local_molecules, molecule_kind_set, molecule_set, region,& 00372 simpar, shell, sum_of_thermostats, para_env, error) 00373 INTEGER, INTENT(IN) :: dis_type 00374 INTEGER, INTENT(OUT) :: natoms_local, nmol_local 00375 INTEGER, DIMENSION(:), POINTER :: const_mol, tot_const 00376 INTEGER, DIMENSION(:, :), POINTER :: point 00377 TYPE(distribution_1d_type), POINTER :: local_molecules 00378 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 00379 TYPE(molecule_type), POINTER :: molecule_set(:) 00380 INTEGER, INTENT(IN) :: region 00381 TYPE(simpar_type), POINTER :: simpar 00382 LOGICAL, INTENT(IN) :: shell 00383 INTEGER, INTENT(IN) :: sum_of_thermostats 00384 TYPE(cp_para_env_type), POINTER :: para_env 00385 TYPE(cp_error_type), INTENT(inout) :: error 00386 00387 CHARACTER(LEN=*), PARAMETER :: routineN = 'adiabatic_region_evaluate', 00388 routineP = moduleN//':'//routineN 00389 00390 INTEGER :: atm_offset, first_atom, handle, icount, ikind, ilist, imol, 00391 imol_local, katom, last_atom, natom, nc, nfixd, nkind, nmol_per_kind, 00392 nmolecule, nshell, stat 00393 LOGICAL :: failure 00394 TYPE(fixd_constraint_type), 00395 DIMENSION(:), POINTER :: fixd_list 00396 TYPE(molecule_kind_type), POINTER :: molecule_kind 00397 TYPE(molecule_type), POINTER :: molecule 00398 00399 CALL timeset(routineN,handle) 00400 00401 natoms_local = 0 00402 nmol_local = 0 00403 nkind = SIZE ( molecule_kind_set ) 00404 NULLIFY(fixd_list, molecule_kind, molecule) 00405 ! Compute the TOTAL number of molecules and atoms on THIS PROC and 00406 ! TOTAL number of molecules of IKIND on THIS PROC 00407 DO ikind = 1, nkind 00408 molecule_kind => molecule_kind_set ( ikind ) 00409 CALL get_molecule_kind ( molecule_kind, natom=natom, nshell=nshell ) 00410 IF (shell) THEN 00411 IF (nshell/=0) THEN 00412 natoms_local = natoms_local + nshell * local_molecules % n_el ( ikind ) 00413 nmol_local = nmol_local + local_molecules % n_el ( ikind ) 00414 END IF 00415 ELSE 00416 natoms_local = natoms_local + natom * local_molecules % n_el ( ikind ) 00417 nmol_local = nmol_local + local_molecules % n_el ( ikind ) 00418 END IF 00419 END DO 00420 00421 CPPostcondition(.NOT.ASSOCIATED(const_mol),cp_failure_level,routineP,error,failure) 00422 CPPostcondition(.NOT.ASSOCIATED(tot_const),cp_failure_level,routineP,error,failure) 00423 CPPostcondition(.NOT.ASSOCIATED(point),cp_failure_level,routineP,error,failure) 00424 IF ( dis_type == do_thermo_no_communication ) THEN 00425 ALLOCATE ( const_mol (nmol_local), STAT = stat ) 00426 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00427 ALLOCATE ( tot_const (nmol_local), STAT = stat ) 00428 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00429 ALLOCATE ( point (2, nmol_local), STAT = stat ) 00430 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00431 00432 point (:,:)= 0 00433 atm_offset = 0 00434 icount = 0 00435 DO ikind = 1, nkind 00436 nmol_per_kind = local_molecules % n_el ( ikind ) 00437 molecule_kind => molecule_kind_set ( ikind ) 00438 CALL get_molecule_kind ( molecule_kind, nconstraint=nc, natom = natom,& 00439 fixd_list=fixd_list, nshell=nshell) 00440 IF (shell) natom = nshell 00441 DO imol_local = 1, nmol_per_kind 00442 icount = icount + 1 00443 point ( 1, icount ) = atm_offset + 1 00444 point ( 2, icount ) = atm_offset + natom 00445 IF (.NOT.shell) THEN 00446 ! nc keeps track of all constraints but not fixed ones.. 00447 ! Let's identify fixed atoms for this molecule 00448 nfixd = 0 00449 imol = local_molecules%list(ikind)%array(imol_local) 00450 molecule => molecule_set(imol) 00451 CALL get_molecule ( molecule, first_atom=first_atom, last_atom=last_atom) 00452 IF (ASSOCIATED(fixd_list)) THEN 00453 DO katom = first_atom, last_atom 00454 DO ilist = 1, SIZE(fixd_list) 00455 IF ( ( katom == fixd_list(ilist)%fixd ) .AND. & 00456 (.NOT. fixd_list(ilist)%restraint%active)) THEN 00457 SELECT CASE(fixd_list(ilist)%itype) 00458 CASE(use_perd_x,use_perd_y,use_perd_z) 00459 nfixd=nfixd+1 00460 CASE(use_perd_xy,use_perd_xz,use_perd_yz) 00461 nfixd=nfixd+2 00462 CASE(use_perd_xyz) 00463 nfixd=nfixd+3 00464 END SELECT 00465 END IF 00466 END DO 00467 END DO 00468 END IF 00469 const_mol ( icount ) = nc + nfixd 00470 tot_const ( icount ) = const_mol ( icount ) 00471 END IF 00472 atm_offset = point ( 2, icount ) 00473 END DO 00474 END DO 00475 ELSE IF ( dis_type == do_thermo_communication ) THEN 00476 ALLOCATE ( const_mol ( nkind ), STAT = stat ) 00477 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00478 ALLOCATE ( tot_const ( nkind ), STAT = stat ) 00479 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00480 ALLOCATE ( point ( 2, nkind ), STAT = stat ) 00481 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00482 point ( :, : ) = 0 00483 atm_offset = 0 00484 ! nc keeps track of all constraints but not fixed ones.. 00485 DO ikind = 1, nkind 00486 nmol_per_kind = local_molecules % n_el ( ikind ) 00487 molecule_kind => molecule_kind_set ( ikind ) 00488 CALL get_molecule_kind ( molecule_kind, nconstraint=nc, natom = natom,& 00489 nmolecule=nmolecule, nconstraint_fixd=nfixd,nshell=nshell) 00490 IF (shell) natom = nshell 00491 IF (.NOT.shell) THEN 00492 const_mol ( ikind ) = nc 00493 ! Let's consider the fixed atoms only for the total number of constraints 00494 ! in case we are in REPLICATED/INTERACTING thermostats 00495 tot_const ( ikind ) = const_mol ( ikind ) * nmolecule + nfixd 00496 END IF 00497 point ( 1, ikind ) = atm_offset + 1 00498 point ( 2, ikind ) = atm_offset + natom * nmol_per_kind 00499 atm_offset = point (2, ikind) 00500 END DO 00501 ENDIF 00502 IF (( .NOT. simpar % constraint ).OR.shell) THEN 00503 const_mol = 0 00504 tot_const = 0 00505 END IF 00506 00507 CALL timestop(handle) 00508 00509 END SUBROUTINE adiabatic_region_evaluate 00510 00511 ! ***************************************************************************** 00515 SUBROUTINE thermostat_mapping_region ( map_info, deg_of_freedom, massive_atom_list,& 00516 molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, simpar, & 00517 number, region, gci, shell, map_loc_thermo_gen, sum_of_thermostats, error) 00518 00519 TYPE(map_info_type), POINTER :: map_info 00520 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, 00521 massive_atom_list 00522 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 00523 TYPE(distribution_1d_type), POINTER :: local_molecules 00524 TYPE(molecule_type), POINTER :: molecule_set(:) 00525 TYPE(cp_para_env_type), POINTER :: para_env 00526 INTEGER, INTENT(OUT) :: natoms_local 00527 TYPE(simpar_type), POINTER :: simpar 00528 INTEGER, INTENT(IN) :: number, region 00529 TYPE(global_constraint_type), POINTER :: gci 00530 LOGICAL, INTENT(IN) :: shell 00531 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen 00532 INTEGER, INTENT(IN) :: sum_of_thermostats 00533 TYPE(cp_error_type), INTENT(inout) :: error 00534 00535 CHARACTER(LEN=*), PARAMETER :: routineN = 'thermostat_mapping_region', 00536 routineP = moduleN//':'//routineN 00537 00538 INTEGER :: handle, nkind, nmol_local, 00539 nsize, number_of_thermostats, 00540 stat 00541 INTEGER, DIMENSION(:), POINTER :: const_mol, tot_const 00542 INTEGER, DIMENSION(:, :), POINTER :: point 00543 LOGICAL :: check, failure 00544 00545 CALL timeset(routineN,handle) 00546 00547 failure = .FALSE. 00548 NULLIFY ( const_mol, tot_const, point) 00549 CPPostcondition(.NOT.ASSOCIATED(deg_of_freedom),cp_failure_level,routineP,error,failure) 00550 CPPostcondition(.NOT.ASSOCIATED(massive_atom_list),cp_failure_level,routineP,error,failure) 00551 00552 nkind = SIZE(molecule_kind_set) 00553 CALL mapping_region_evaluate(map_info%dis_type, natoms_local, nmol_local,& 00554 const_mol, tot_const, point, local_molecules, molecule_kind_set, molecule_set,& 00555 region, simpar, shell, map_loc_thermo_gen, sum_of_thermostats, para_env, error) 00556 00557 ! Now we can allocate the target array s_kin and p_kin.. 00558 SELECT CASE(region) 00559 CASE(do_region_global, do_region_molecule, do_region_massive) 00560 nsize = number 00561 CASE(do_region_defined) 00562 nsize = sum_of_thermostats 00563 END SELECT 00564 ALLOCATE (map_info%s_kin(nsize),STAT=stat) 00565 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00566 ALLOCATE (map_info%v_scale(nsize),STAT=stat) 00567 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00568 ALLOCATE (map_info%p_kin(3,natoms_local),STAT=stat) 00569 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00570 ALLOCATE (map_info%p_scale(3,natoms_local),STAT=stat) 00571 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00572 ! Allocate index array 00573 ALLOCATE ( map_info%index(number), STAT=stat) 00574 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00575 ALLOCATE ( map_info%map_index(number), STAT=stat) 00576 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00577 ALLOCATE (deg_of_freedom(number),STAT=stat) 00578 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00579 00580 CALL massive_list_generate (molecule_set, molecule_kind_set, & 00581 local_molecules, para_env, massive_atom_list, region, shell,& 00582 error) 00583 00584 CALL thermostat_mapping_region_low(region, map_info, nkind, point,& 00585 deg_of_freedom, local_molecules, const_mol, massive_atom_list,& 00586 tot_const, molecule_set, number_of_thermostats, shell, gci,& 00587 map_loc_thermo_gen, error) 00588 00589 check = (number==number_of_thermostats) 00590 CPPrecondition(check,cp_fatal_level,routineP,error,failure) 00591 DEALLOCATE (const_mol,STAT=stat) 00592 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00593 DEALLOCATE (tot_const,STAT=stat) 00594 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00595 DEALLOCATE (point, STAT = stat ) 00596 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00597 00598 CALL timestop(handle) 00599 00600 END SUBROUTINE thermostat_mapping_region 00601 00602 ! ***************************************************************************** 00606 SUBROUTINE thermostat_mapping_region_low(region, map_info, nkind, point,& 00607 deg_of_freedom, local_molecules, const_mol, massive_atom_list, tot_const,& 00608 molecule_set, number, shell, gci, map_loc_thermo_gen,error) 00609 00610 INTEGER, INTENT(IN) :: region 00611 TYPE(map_info_type), POINTER :: map_info 00612 INTEGER :: nkind 00613 INTEGER, DIMENSION(:, :), POINTER :: point 00614 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom 00615 TYPE(distribution_1d_type), POINTER :: local_molecules 00616 INTEGER, DIMENSION(:), POINTER :: const_mol, massive_atom_list, 00617 tot_const 00618 TYPE(molecule_type), POINTER :: molecule_set(:) 00619 INTEGER, INTENT(OUT) :: number 00620 LOGICAL, INTENT(IN) :: shell 00621 TYPE(global_constraint_type), POINTER :: gci 00622 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen 00623 TYPE(cp_error_type), INTENT(inout) :: error 00624 00625 CHARACTER(LEN=*), PARAMETER :: 00626 routineN = 'thermostat_mapping_region_low', 00627 routineP = moduleN//':'//routineN 00628 00629 INTEGER :: first_atom, first_shell, handle, i, icount, ielement, ii, 00630 ikind, imap, imol, imol_local, ipart, itmp, jj, k, kk, last_atom, 00631 last_shell, nglob_cns, nmol_local, stat 00632 INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp, wrk 00633 LOGICAL :: check, failure, 00634 global_constraints 00635 TYPE(molecule_type), POINTER :: molecule 00636 00637 CALL timeset(routineN,handle) 00638 00639 failure = .FALSE. 00640 global_constraints = ASSOCIATED(gci) 00641 deg_of_freedom = 0 00642 icount = 0 00643 number = 0 00644 nglob_cns = 0 00645 IF (global_constraints) nglob_cns = gci%ntot-gci%nrestraint 00646 IF ( region==do_region_global) THEN 00647 ! Global Region 00648 check = ( map_info%dis_type == do_thermo_communication ) 00649 CPPostcondition(check,cp_failure_level,routineP,error,failure) 00650 DO ikind = 1, nkind 00651 DO jj = point ( 1, ikind ), point ( 2, ikind ) 00652 DO ii = 1, 3 00653 map_info%p_kin(ii,jj)%point => map_info%s_kin(1) 00654 map_info%p_scale(ii,jj)%point => map_info%v_scale(1) 00655 END DO 00656 END DO 00657 deg_of_freedom(1) = deg_of_freedom(1) + tot_const(ikind) 00658 map_info%index(1) = 1 00659 map_info%map_index(1) = 1 00660 number = 1 00661 END DO 00662 deg_of_freedom(1) = deg_of_freedom(1) + nglob_cns 00663 ELSE IF ( region==do_region_defined) THEN 00664 ! User defined Region to thermostat 00665 check = ( map_info%dis_type == do_thermo_communication ) 00666 CPPostcondition(check,cp_failure_level,routineP,error,failure) 00667 ! Lets' identify the matching of the local thermostat w.r.t. the global one 00668 itmp = SIZE(map_loc_thermo_gen) 00669 ALLOCATE(tmp(itmp),stat=stat) 00670 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00671 ALLOCATE(wrk(itmp),stat=stat) 00672 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00673 tmp = map_loc_thermo_gen 00674 CALL sort(tmp, itmp, wrk) 00675 number = 1 00676 map_info%index(number) = tmp(1) 00677 map_info%map_index(number) = tmp(1) 00678 deg_of_freedom(number) = tot_const(tmp(1)) 00679 DO i = 2, itmp 00680 IF (tmp(i)/=tmp(i-1)) THEN 00681 number = number + 1 00682 map_info%index(number) = tmp(i) 00683 map_info%map_index(number) = tmp(i) 00684 deg_of_freedom(number) = tot_const(tmp(i)) 00685 END IF 00686 END DO 00687 DEALLOCATE(tmp,stat=stat) 00688 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00689 DEALLOCATE(wrk,stat=stat) 00690 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00691 DO jj = 1, SIZE(map_loc_thermo_gen) 00692 DO ii = 1, 3 00693 imap = map_loc_thermo_gen(jj) 00694 map_info%p_kin(ii,jj)%point => map_info%s_kin(imap) 00695 map_info%p_scale(ii,jj)%point => map_info%v_scale(imap) 00696 END DO 00697 END DO 00698 IF (nglob_cns/=0) THEN 00699 CALL stop_program(routineN,moduleN,__LINE__,& 00700 "User Defined thermostats with global constraints not implemented!") 00701 END IF 00702 ELSE IF ( region==do_region_molecule) THEN 00703 ! Molecular Region 00704 IF ( map_info%dis_type == do_thermo_no_communication ) THEN 00705 ! This is the standard case.. 00706 DO ikind = 1, nkind 00707 nmol_local = local_molecules % n_el (ikind) 00708 DO imol_local = 1, nmol_local 00709 imol = local_molecules%list(ikind)%array(imol_local) 00710 number = number + 1 00711 map_info%index(number) = imol 00712 map_info%map_index(number) = number 00713 deg_of_freedom ( number ) = const_mol ( number ) 00714 DO kk = point ( 1, number ), point ( 2, number ) 00715 DO jj = 1, 3 00716 map_info%p_kin(jj,kk) %point => map_info%s_kin(number) 00717 map_info%p_scale(jj,kk) %point => map_info%v_scale(number) 00718 END DO 00719 END DO 00720 END DO 00721 END DO 00722 ELSE IF ( map_info%dis_type == do_thermo_communication ) THEN 00723 ! This case is quite rare and happens only when we have one molecular 00724 ! kind and one molecule.. 00725 CPPostcondition(nkind==1,cp_failure_level,routineP,error,failure) 00726 number = number + 1 00727 map_info%index(number) = number 00728 map_info%map_index(number) = number 00729 deg_of_freedom ( number ) = deg_of_freedom ( number ) + tot_const ( nkind ) 00730 DO kk = point ( 1, nkind ), point ( 2, nkind ) 00731 DO jj = 1, 3 00732 map_info%p_kin ( jj, kk ) % point => map_info%s_kin(number) 00733 map_info%p_scale ( jj, kk ) % point => map_info%v_scale(number) 00734 END DO 00735 END DO 00736 ELSE 00737 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 00738 END IF 00739 IF (nglob_cns/=0) THEN 00740 CALL stop_program(routineN,moduleN,__LINE__,& 00741 "Molecular thermostats with global constraints are impossible!") 00742 END IF 00743 ELSE IF ( region==do_region_massive) THEN 00744 ! Massive Region 00745 check = ( map_info%dis_type == do_thermo_no_communication ) 00746 CPPostcondition(check,cp_failure_level,routineP,error,failure) 00747 DO ikind = 1, nkind 00748 nmol_local = local_molecules % n_el ( ikind ) 00749 DO imol_local = 1, nmol_local 00750 icount = icount + 1 00751 imol = local_molecules % list ( ikind ) % array ( imol_local ) 00752 molecule => molecule_set ( imol ) 00753 CALL get_molecule ( molecule, first_atom=first_atom, last_atom=last_atom,& 00754 first_shell=first_shell, last_shell=last_shell) 00755 IF (shell) THEN 00756 first_atom = first_shell 00757 last_atom = last_shell 00758 ELSE 00759 IF ((tot_const(icount)>0).OR.(nglob_cns/=0)) THEN 00760 CALL stop_program(routineN,moduleN,__LINE__,& 00761 "Massive thermostats with constraints are impossible!") 00762 END IF 00763 END IF 00764 k = 0 00765 DO ii = point ( 1, icount ), point ( 2, icount ) 00766 ipart = first_atom + k 00767 ielement = locate(massive_atom_list,ipart) 00768 k = k + 1 00769 DO jj = 1, 3 00770 number = number + 1 00771 map_info%index(number) = (ielement - 1)*3 + jj 00772 map_info%map_index(number) = number 00773 map_info%p_kin(jj,ii) %point => map_info%s_kin(number) 00774 map_info%p_scale(jj,ii) %point => map_info%v_scale(number) 00775 END DO 00776 END DO 00777 IF ( first_atom + k -1 /= last_atom ) THEN 00778 CALL stop_program(routineN,moduleN,__LINE__,& 00779 "Inconsistent mapping of particles") 00780 END IF 00781 END DO 00782 END DO 00783 ELSE 00784 CALL stop_program(routineN,moduleN,__LINE__,"Invalid region!") 00785 END IF 00786 00787 CALL timestop(handle) 00788 00789 END SUBROUTINE thermostat_mapping_region_low 00790 00791 ! ***************************************************************************** 00795 SUBROUTINE mapping_region_evaluate(dis_type, natoms_local, nmol_local, const_mol,& 00796 tot_const, point, local_molecules, molecule_kind_set, molecule_set, region,& 00797 simpar, shell, map_loc_thermo_gen, sum_of_thermostats, para_env, error) 00798 INTEGER, INTENT(IN) :: dis_type 00799 INTEGER, INTENT(OUT) :: natoms_local, nmol_local 00800 INTEGER, DIMENSION(:), POINTER :: const_mol, tot_const 00801 INTEGER, DIMENSION(:, :), POINTER :: point 00802 TYPE(distribution_1d_type), POINTER :: local_molecules 00803 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 00804 TYPE(molecule_type), POINTER :: molecule_set(:) 00805 INTEGER, INTENT(IN) :: region 00806 TYPE(simpar_type), POINTER :: simpar 00807 LOGICAL, INTENT(IN) :: shell 00808 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen 00809 INTEGER, INTENT(IN) :: sum_of_thermostats 00810 TYPE(cp_para_env_type), POINTER :: para_env 00811 TYPE(cp_error_type), INTENT(inout) :: error 00812 00813 CHARACTER(LEN=*), PARAMETER :: routineN = 'mapping_region_evaluate', 00814 routineP = moduleN//':'//routineN 00815 00816 INTEGER :: atm_offset, first_atom, handle, i, iatm, icount, id_region, 00817 ikind, ilist, imol, imol_local, j, jatm, katom, last_atom, natom, nc, 00818 nfixd, nkind, nmol_per_kind, nmolecule, nshell, stat 00819 LOGICAL :: failure 00820 TYPE(colvar_constraint_type), 00821 DIMENSION(:), POINTER :: colv_list 00822 TYPE(fixd_constraint_type), 00823 DIMENSION(:), POINTER :: fixd_list 00824 TYPE(g3x3_constraint_type), 00825 DIMENSION(:), POINTER :: g3x3_list 00826 TYPE(g4x6_constraint_type), 00827 DIMENSION(:), POINTER :: g4x6_list 00828 TYPE(molecule_kind_type), POINTER :: molecule_kind 00829 TYPE(molecule_type), POINTER :: molecule 00830 00831 CALL timeset(routineN,handle) 00832 00833 natoms_local = 0 00834 nmol_local = 0 00835 nkind = SIZE ( molecule_kind_set ) 00836 NULLIFY(fixd_list, molecule_kind, molecule, colv_list, g3x3_list, g4x6_list) 00837 ! Compute the TOTAL number of molecules and atoms on THIS PROC and 00838 ! TOTAL number of molecules of IKIND on THIS PROC 00839 DO ikind = 1, nkind 00840 molecule_kind => molecule_kind_set ( ikind ) 00841 CALL get_molecule_kind ( molecule_kind, natom=natom, nshell=nshell ) 00842 IF (shell) THEN 00843 IF (nshell/=0) THEN 00844 natoms_local = natoms_local + nshell * local_molecules % n_el ( ikind ) 00845 nmol_local = nmol_local + local_molecules % n_el ( ikind ) 00846 END IF 00847 ELSE 00848 natoms_local = natoms_local + natom * local_molecules % n_el ( ikind ) 00849 nmol_local = nmol_local + local_molecules % n_el ( ikind ) 00850 END IF 00851 END DO 00852 00853 CPPostcondition(.NOT.ASSOCIATED(const_mol),cp_failure_level,routineP,error,failure) 00854 CPPostcondition(.NOT.ASSOCIATED(tot_const),cp_failure_level,routineP,error,failure) 00855 CPPostcondition(.NOT.ASSOCIATED(point),cp_failure_level,routineP,error,failure) 00856 IF ( dis_type == do_thermo_no_communication ) THEN 00857 ALLOCATE ( const_mol (nmol_local), STAT = stat ) 00858 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00859 ALLOCATE ( tot_const (nmol_local), STAT = stat ) 00860 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00861 ALLOCATE ( point (2, nmol_local), STAT = stat ) 00862 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00863 00864 point (:,:)= 0 00865 atm_offset = 0 00866 icount = 0 00867 DO ikind = 1, nkind 00868 nmol_per_kind = local_molecules % n_el ( ikind ) 00869 molecule_kind => molecule_kind_set ( ikind ) 00870 CALL get_molecule_kind ( molecule_kind, nconstraint=nc, natom = natom,& 00871 fixd_list=fixd_list, nshell=nshell) 00872 IF (shell) natom = nshell 00873 DO imol_local = 1, nmol_per_kind 00874 icount = icount + 1 00875 point ( 1, icount ) = atm_offset + 1 00876 point ( 2, icount ) = atm_offset + natom 00877 IF (.NOT.shell) THEN 00878 ! nc keeps track of all constraints but not fixed ones.. 00879 ! Let's identify fixed atoms for this molecule 00880 nfixd = 0 00881 imol = local_molecules%list(ikind)%array(imol_local) 00882 molecule => molecule_set(imol) 00883 CALL get_molecule ( molecule, first_atom=first_atom, last_atom=last_atom) 00884 IF (ASSOCIATED(fixd_list)) THEN 00885 DO katom = first_atom, last_atom 00886 DO ilist = 1, SIZE(fixd_list) 00887 IF ( ( katom == fixd_list(ilist)%fixd ) .AND. & 00888 (.NOT. fixd_list(ilist)%restraint%active)) THEN 00889 SELECT CASE(fixd_list(ilist)%itype) 00890 CASE(use_perd_x,use_perd_y,use_perd_z) 00891 nfixd=nfixd+1 00892 CASE(use_perd_xy,use_perd_xz,use_perd_yz) 00893 nfixd=nfixd+2 00894 CASE(use_perd_xyz) 00895 nfixd=nfixd+3 00896 END SELECT 00897 END IF 00898 END DO 00899 END DO 00900 END IF 00901 const_mol ( icount ) = nc + nfixd 00902 tot_const ( icount ) = const_mol ( icount ) 00903 END IF 00904 atm_offset = point ( 2, icount ) 00905 END DO 00906 END DO 00907 ELSE IF ( dis_type == do_thermo_communication ) THEN 00908 IF (region==do_region_defined) THEN 00909 ! Setup of the arbitrary region 00910 ALLOCATE ( tot_const (sum_of_thermostats), STAT = stat ) 00911 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00912 ALLOCATE ( point (2, 0), STAT = stat ) 00913 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00914 ALLOCATE ( const_mol (0), STAT = stat ) 00915 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00916 atm_offset = 0 00917 tot_const = 0 00918 const_mol = 0 00919 point = 0 00920 DO ikind = 1, nkind 00921 nmol_per_kind = local_molecules % n_el ( ikind ) 00922 molecule_kind => molecule_kind_set ( ikind ) 00923 CALL get_molecule_kind ( molecule_kind, nconstraint=nc, natom = natom,& 00924 fixd_list=fixd_list, colv_list=colv_list, g3x3_list=g3x3_list,& 00925 g4x6_list=g4x6_list, nshell=nshell) 00926 IF (shell) natom = nshell 00927 DO imol_local = 1, nmol_per_kind 00928 IF (.NOT.shell) THEN 00929 ! First if nc is not zero let's check if all atoms of a molecule 00930 ! are in the same thermostatting region.. 00931 imol = local_molecules%list(ikind)%array(imol_local) 00932 molecule => molecule_set(imol) 00933 id_region = map_loc_thermo_gen(atm_offset+1) 00934 IF (ALL(map_loc_thermo_gen(atm_offset+1:atm_offset+natom)==id_region)) THEN 00935 ! All the atoms of a molecule are within the same thermostatting 00936 ! region.. this is the easy case.. 00937 tot_const(id_region) = tot_const(id_region) + nc 00938 ELSE 00939 ! If not let's check the single constraints defined for this molecule 00940 ! and continue only when atoms involved in the constraint belong to 00941 ! the same thermostatting region 00942 IF (ASSOCIATED(colv_list)) THEN 00943 DO i = 1, SIZE(colv_list) 00944 IF (.NOT.colv_list(i)%restraint%active) THEN 00945 iatm = atm_offset + colv_list(i)%i_atoms(1) 00946 DO j = 2, SIZE(colv_list(i)%i_atoms) 00947 jatm = atm_offset + colv_list(i)%i_atoms(j) 00948 IF (map_loc_thermo_gen(iatm)/=map_loc_thermo_gen(jatm)) THEN 00949 CALL stop_program(routineN,moduleN,__LINE__,& 00950 "User Defined Region: "//& 00951 "A constraint (COLV) was defined between two thermostatting regions! "//& 00952 "This is not allowed!") 00953 END IF 00954 END DO 00955 id_region = map_loc_thermo_gen(iatm) 00956 tot_const(id_region) = tot_const(id_region) + 1 00957 END IF 00958 END DO 00959 END IF 00960 IF (ASSOCIATED(g3x3_list)) THEN 00961 DO i = 1, SIZE(g3x3_list) 00962 IF (.NOT.g3x3_list(i)%restraint%active) THEN 00963 iatm = atm_offset + g3x3_list(i)%a 00964 jatm = atm_offset + g3x3_list(i)%b 00965 IF (map_loc_thermo_gen(iatm)/=map_loc_thermo_gen(jatm)) THEN 00966 CALL stop_program(routineN,moduleN,__LINE__,& 00967 "User Defined Region: "//& 00968 "A constraint (G3X3) was defined between two thermostatting regions! "//& 00969 "This is not allowed!") 00970 END IF 00971 jatm = atm_offset + g3x3_list(i)%c 00972 IF (map_loc_thermo_gen(iatm)/=map_loc_thermo_gen(jatm)) THEN 00973 CALL stop_program(routineN,moduleN,__LINE__,& 00974 "User Defined Region: "//& 00975 "A constraint (G3X3) was defined between two thermostatting regions! "//& 00976 "This is not allowed!") 00977 END IF 00978 END IF 00979 id_region = map_loc_thermo_gen(iatm) 00980 tot_const(id_region) = tot_const(id_region) + 3 00981 END DO 00982 END IF 00983 IF (ASSOCIATED(g4x6_list)) THEN 00984 DO i = 1, SIZE(g4x6_list) 00985 IF (.NOT.g4x6_list(i)%restraint%active) THEN 00986 iatm = atm_offset + g4x6_list(i)%a 00987 jatm = atm_offset + g4x6_list(i)%b 00988 IF (map_loc_thermo_gen(iatm)/=map_loc_thermo_gen(jatm)) THEN 00989 CALL stop_program(routineN,moduleN,__LINE__,& 00990 " User Defined Region: "//& 00991 "A constraint (G4X6) was defined between two thermostatting regions! "//& 00992 "This is not allowed!") 00993 END IF 00994 jatm = atm_offset + g4x6_list(i)%c 00995 IF (map_loc_thermo_gen(iatm)/=map_loc_thermo_gen(jatm)) THEN 00996 CALL stop_program(routineN,moduleN,__LINE__,& 00997 " User Defined Region: "//& 00998 "A constraint (G4X6) was defined between two thermostatting regions! "//& 00999 "This is not allowed!") 01000 END IF 01001 jatm = atm_offset + g4x6_list(i)%d 01002 IF (map_loc_thermo_gen(iatm)/=map_loc_thermo_gen(jatm)) THEN 01003 CALL stop_program(routineN,moduleN,__LINE__,& 01004 " User Defined Region: "//& 01005 "A constraint (G4X6) was defined between two thermostatting regions! "//& 01006 "This is not allowed!") 01007 END IF 01008 END IF 01009 id_region = map_loc_thermo_gen(iatm) 01010 tot_const(id_region) = tot_const(id_region) + 6 01011 END DO 01012 END IF 01013 END IF 01014 ! Here we handle possibly fixed atoms 01015 IF (ASSOCIATED(fixd_list)) THEN 01016 CALL get_molecule ( molecule, first_atom=first_atom, last_atom=last_atom) 01017 iatm = 0 01018 DO katom = first_atom, last_atom 01019 iatm = iatm + 1 01020 DO ilist = 1, SIZE(fixd_list) 01021 IF ( ( katom == fixd_list(ilist)%fixd ) .AND. & 01022 (.NOT. fixd_list(ilist)%restraint%active)) THEN 01023 id_region = map_loc_thermo_gen(atm_offset+iatm) 01024 SELECT CASE(fixd_list(ilist)%itype) 01025 CASE(use_perd_x, use_perd_y, use_perd_z) 01026 tot_const(id_region) = tot_const(id_region) + 1 01027 CASE(use_perd_xy, use_perd_xz, use_perd_yz) 01028 tot_const(id_region) = tot_const(id_region) + 2 01029 CASE(use_perd_xyz) 01030 tot_const(id_region) = tot_const(id_region) + 3 01031 END SELECT 01032 END IF 01033 END DO 01034 END DO 01035 END IF 01036 END IF 01037 atm_offset = atm_offset + natom 01038 END DO 01039 END DO 01040 CALL mp_sum(tot_const, para_env%group) 01041 ELSE 01042 ALLOCATE ( const_mol ( nkind ), STAT = stat ) 01043 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01044 ALLOCATE ( tot_const ( nkind ), STAT = stat ) 01045 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01046 ALLOCATE ( point ( 2, nkind ), STAT = stat ) 01047 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01048 point ( :, : ) = 0 01049 atm_offset = 0 01050 ! nc keeps track of all constraints but not fixed ones.. 01051 DO ikind = 1, nkind 01052 nmol_per_kind = local_molecules % n_el ( ikind ) 01053 molecule_kind => molecule_kind_set ( ikind ) 01054 CALL get_molecule_kind ( molecule_kind, nconstraint=nc, natom = natom,& 01055 nmolecule=nmolecule, nconstraint_fixd=nfixd,nshell=nshell) 01056 IF (shell) natom = nshell 01057 IF (.NOT.shell) THEN 01058 const_mol ( ikind ) = nc 01059 ! Let's consider the fixed atoms only for the total number of constraints 01060 ! in case we are in REPLICATED/INTERACTING thermostats 01061 tot_const ( ikind ) = const_mol ( ikind ) * nmolecule + nfixd 01062 END IF 01063 point ( 1, ikind ) = atm_offset + 1 01064 point ( 2, ikind ) = atm_offset + natom * nmol_per_kind 01065 atm_offset = point (2, ikind) 01066 END DO 01067 END IF 01068 ENDIF 01069 IF (( .NOT. simpar % constraint ).OR.shell) THEN 01070 const_mol = 0 01071 tot_const = 0 01072 END IF 01073 01074 CALL timestop(handle) 01075 01076 END SUBROUTINE mapping_region_evaluate 01077 01078 ! ***************************************************************************** 01079 SUBROUTINE massive_list_generate ( molecule_set, molecule_kind_set, & 01080 local_molecules, para_env, massive_atom_list, region, shell, error ) 01081 01082 TYPE(molecule_type), POINTER :: molecule_set( : ) 01083 TYPE(molecule_kind_type), POINTER :: molecule_kind_set( : ) 01084 TYPE(distribution_1d_type), POINTER :: local_molecules 01085 TYPE(cp_para_env_type), POINTER :: para_env 01086 INTEGER, POINTER :: massive_atom_list( : ) 01087 INTEGER, INTENT(IN) :: region 01088 LOGICAL, INTENT(IN) :: shell 01089 TYPE(cp_error_type), INTENT(inout) :: error 01090 01091 CHARACTER(LEN=*), PARAMETER :: routineN = 'massive_list_generate', 01092 routineP = moduleN//':'//routineN 01093 01094 INTEGER :: first_atom, first_shell, handle, i, ikind, imol, iproc, j, 01095 natom, ncount, nkind, nmol_per_kind, nshell, num_massive_atm, 01096 num_massive_atm_local, offset, stat 01097 INTEGER, DIMENSION(:), POINTER :: array_num_massive_atm, 01098 local_atm_list, work 01099 LOGICAL :: failure 01100 TYPE(molecule_kind_type), POINTER :: molecule_kind 01101 TYPE(molecule_type), POINTER :: molecule 01102 01103 CALL timeset(routineN,handle) 01104 01105 failure = .FALSE. 01106 num_massive_atm_local = 0 01107 NULLIFY(local_atm_list) 01108 CALL reallocate(local_atm_list,1,num_massive_atm_local) 01109 01110 nkind = SIZE ( molecule_kind_set ) 01111 DO ikind = 1, nkind 01112 nmol_per_kind = local_molecules%n_el(ikind) 01113 DO imol = 1, nmol_per_kind 01114 i = local_molecules%list(ikind)%array(imol) 01115 molecule => molecule_set ( i ) 01116 molecule_kind => molecule % molecule_kind 01117 CALL get_molecule_kind(molecule_kind,natom=natom,nshell=nshell) 01118 IF(region == do_region_massive) THEN 01119 IF (shell) THEN 01120 natom = nshell 01121 END IF 01122 num_massive_atm_local = num_massive_atm_local + natom 01123 CALL reallocate(local_atm_list,1,num_massive_atm_local) 01124 CALL get_molecule (molecule,first_atom=first_atom,first_shell=first_shell) 01125 IF (shell) THEN 01126 first_atom = first_shell 01127 END IF 01128 DO j=1,natom 01129 local_atm_list(num_massive_atm_local-natom+j) = first_atom -1 + j 01130 END DO 01131 END IF 01132 END DO 01133 END DO 01134 01135 ALLOCATE(array_num_massive_atm(para_env%num_pe), STAT = stat ) 01136 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01137 CALL mp_allgather(num_massive_atm_local,array_num_massive_atm,para_env%group) 01138 01139 num_massive_atm = SUM(array_num_massive_atm) 01140 ALLOCATE(massive_atom_list(num_massive_atm), STAT = stat ) 01141 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01142 01143 offset = 0 01144 DO iproc=1,para_env%num_pe 01145 ncount = array_num_massive_atm(iproc) 01146 ALLOCATE(work(ncount), STAT = stat ) 01147 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01148 IF(para_env%mepos == (iproc-1)) THEN 01149 DO i=1,ncount 01150 work(i) = local_atm_list(i) 01151 END DO 01152 ELSE 01153 work(:) = 0 01154 END IF 01155 CALL mp_bcast(work,iproc-1,para_env%group) 01156 DO i=1,ncount 01157 massive_atom_list(offset+i) = work(i) 01158 END DO 01159 DEALLOCATE(work, STAT = stat) 01160 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01161 offset = offset + array_num_massive_atm(iproc) 01162 END DO 01163 01164 ! Sort atom list 01165 ALLOCATE (work(num_massive_atm),STAT=stat) 01166 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01167 "work",int_size*num_massive_atm) 01168 CALL sort(massive_atom_list,num_massive_atm,work) 01169 DEALLOCATE (work,STAT=stat) 01170 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"work") 01171 01172 DEALLOCATE(local_atm_list, STAT = stat) 01173 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01174 DEALLOCATE(array_num_massive_atm, STAT = stat ) 01175 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01176 01177 CALL timestop(handle) 01178 01179 END SUBROUTINE massive_list_generate 01180 01181 ! ***************************************************************************** 01185 SUBROUTINE init_baro_map_info(map_info, ndeg, num_thermo, error) 01186 01187 TYPE(map_info_type), POINTER :: map_info 01188 INTEGER, INTENT(IN) :: ndeg, num_thermo 01189 TYPE(cp_error_type), INTENT(inout) :: error 01190 01191 CHARACTER(LEN=*), PARAMETER :: routineN = 'init_baro_map_info', 01192 routineP = moduleN//':'//routineN 01193 01194 INTEGER :: handle, i, stat 01195 LOGICAL :: failure 01196 01197 CALL timeset(routineN,handle) 01198 01199 ALLOCATE (map_info%s_kin(num_thermo),stat=stat) 01200 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01201 ALLOCATE (map_info%v_scale(num_thermo),stat=stat ) 01202 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01203 ALLOCATE (map_info%p_kin(1,ndeg),stat=stat ) 01204 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01205 ALLOCATE (map_info%p_scale(1,ndeg),stat=stat) 01206 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01207 ! Allocate the index array 01208 ALLOCATE (map_info%index(1),stat=stat) 01209 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01210 ALLOCATE (map_info%map_index(1),stat=stat) 01211 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01212 01213 ! Begin the mapping loop 01214 DO i = 1, ndeg 01215 map_info%p_kin(1,i)%point => map_info%s_kin(1) 01216 map_info%p_scale(1,i)%point => map_info%v_scale(1) 01217 END DO 01218 map_info%index(1) = 1 01219 map_info%map_index(1) = 1 01220 01221 CALL timestop(handle) 01222 01223 END SUBROUTINE init_baro_map_info 01224 01225 END MODULE thermostat_mapping 01226
1.7.3