CP2K 2.4 (Revision 12889)

wiener_process.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 wiener_process
00014 
00015   USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
00016   USE cp_para_types,                   ONLY: cp_para_env_type
00017   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
00018                                              cp_subsys_type
00019   USE distribution_1d_types,           ONLY: distribution_1d_type,&
00020                                              init_local_particle_set
00021   USE f77_blas
00022   USE force_env_types,                 ONLY: force_env_get,&
00023                                              force_env_type
00024   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
00025                                              section_vals_type
00026   USE kinds,                           ONLY: dp
00027   USE md_environment_types,            ONLY: get_md_env,&
00028                                              md_environment_type,&
00029                                              need_per_atom_wiener_process
00030   USE metadynamics_types,              ONLY: meta_env_type
00031   USE parallel_rng_types,              ONLY: GAUSSIAN,&
00032                                              create_rng_stream,&
00033                                              next_rng_seed
00034   USE particle_list_types,             ONLY: particle_list_type
00035   USE simpar_types,                    ONLY: simpar_type
00036   USE string_utilities,                ONLY: compress
00037 #include "cp_common_uses.h"
00038 
00039   IMPLICIT NONE
00040 
00041   PRIVATE
00042 
00043   ! Global parameters in this module
00044   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'wiener_process'
00045 
00046   ! Public subroutines
00047   PUBLIC :: create_wiener_process, create_wiener_process_cv
00048 
00049 CONTAINS
00050 
00051 ! *****************************************************************************
00058   SUBROUTINE create_wiener_process(md_env,error)
00059 
00060     TYPE(md_environment_type), POINTER       :: md_env
00061     TYPE(cp_error_type), INTENT(INOUT)       :: error
00062 
00063     CHARACTER(LEN=*), PARAMETER :: routineN = 'create_wiener_process', 
00064       routineP = moduleN//':'//routineN
00065 
00066     CHARACTER(LEN=40)                        :: name
00067     INTEGER :: iparticle, iparticle_kind, iparticle_local, nparticle, 
00068       nparticle_kind, nparticle_local, stat
00069     LOGICAL                                  :: failure
00070     REAL(KIND=dp), ALLOCATABLE, 
00071       DIMENSION(:, :, :)                     :: seed
00072     REAL(KIND=dp), DIMENSION(3, 2)           :: initial_seed
00073     TYPE(atomic_kind_list_type), POINTER     :: atomic_kinds
00074     TYPE(cp_para_env_type), POINTER          :: para_env
00075     TYPE(cp_subsys_type), POINTER            :: subsys
00076     TYPE(distribution_1d_type), POINTER      :: local_particles
00077     TYPE(force_env_type), POINTER            :: force_env
00078     TYPE(particle_list_type), POINTER        :: particles
00079     TYPE(section_vals_type), POINTER         :: force_env_section, 
00080                                                 subsys_section, work_section
00081     TYPE(simpar_type), POINTER               :: simpar
00082 
00083     failure = .FALSE.
00084     NULLIFY(work_section,force_env)
00085     CPPrecondition (ASSOCIATED(md_env),cp_failure_level,routineP,error,failure)
00086 
00087     CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env,&
00088          simpar=simpar, error=error)
00089 
00090     ![NB] shouldn't the calling process know if it's needed
00091     IF (need_per_atom_wiener_process(md_env, error=error)) THEN
00092 
00093        ! Load initial seed (not needed for a restart)
00094 
00095        initial_seed = next_rng_seed(error=error)
00096 
00097        CALL force_env_get(force_env,force_env_section=force_env_section,&
00098             subsys=subsys,error=error)
00099 
00100        subsys_section => section_vals_get_subs_vals(force_env_section,"SUBSYS",error=error)
00101 
00102        CALL cp_subsys_get(subsys=subsys,atomic_kinds=atomic_kinds,local_particles=local_particles,&
00103             particles=particles,error=error)
00104 
00105        nparticle_kind = atomic_kinds%n_els
00106        nparticle = particles%n_els
00107 
00108        ! Allocate the (local) data structures for the Wiener process
00109 
00110        ALLOCATE(local_particles%local_particle_set(nparticle_kind),STAT=stat)
00111        CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00112 
00113        DO iparticle_kind=1,nparticle_kind
00114           nparticle_local = local_particles%n_el(iparticle_kind)
00115           ALLOCATE (local_particles%local_particle_set(iparticle_kind)%rng(nparticle_local),STAT=stat)
00116           CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00117           DO iparticle_local=1,nparticle_local
00118              NULLIFY (local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream)
00119           END DO
00120        END DO
00121 
00122        ! Each process generates all seeds. The seed generation should be
00123        ! quite fast and in this way a broadcast is avoided.
00124 
00125        ALLOCATE (seed(3,2,nparticle),STAT=stat)
00126        CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00127 
00128        seed(:,:,1) = initial_seed
00129        DO iparticle=2,nparticle
00130           seed(:,:,iparticle) = next_rng_seed(seed(:,:,iparticle-1),error=error)
00131        END DO
00132 
00133        ! Update initial seed
00134 
00135        initial_seed = next_rng_seed(seed(:,:,nparticle),error=error)
00136 
00137        ! Create a random number stream (Wiener process) for each particle
00138 
00139        DO iparticle_kind=1,nparticle_kind
00140           nparticle_local = local_particles%n_el(iparticle_kind)
00141           DO iparticle_local=1,nparticle_local
00142              iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
00143              WRITE (UNIT=name,FMT="(A,I8)") "Wiener process for particle",iparticle
00144              CALL compress(name)
00145              CALL create_rng_stream(rng_stream=local_particles%local_particle_set(iparticle_kind)%&
00146                   &rng(iparticle_local)%stream,name=name,distribution_type=GAUSSIAN,&
00147                   extended_precision=.TRUE., seed=seed(:,:,iparticle),error=error)
00148           END DO
00149        END DO
00150 
00151        DEALLOCATE (seed,STAT=stat)
00152        CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00153 
00154        ! Possibly restart Wiener process
00155        NULLIFY (work_section)
00156        work_section => section_vals_get_subs_vals(section_vals=subsys_section,&
00157             subsection_name="RNG_INIT", error=error)
00158        CALL init_local_particle_set(distribution_1d=local_particles,&
00159             nparticle_kind=nparticle_kind, &
00160             work_section=work_section,error=error)
00161     END IF
00162 
00163   END SUBROUTINE create_wiener_process
00164 
00165 ! *****************************************************************************
00173   SUBROUTINE create_wiener_process_cv(meta_env,error)
00174 
00175     TYPE(meta_env_type), POINTER             :: meta_env
00176     TYPE(cp_error_type), INTENT(INOUT)       :: error
00177 
00178     CHARACTER(LEN=*), PARAMETER :: routineN = 'create_wiener_process_cv', 
00179       routineP = moduleN//':'//routineN
00180 
00181     CHARACTER(LEN=40)                        :: name
00182     INTEGER                                  :: i_c, stat
00183     LOGICAL                                  :: failure
00184     REAL(KIND=dp), ALLOCATABLE, 
00185       DIMENSION(:, :, :)                     :: seed
00186     REAL(KIND=dp), DIMENSION(3, 2)           :: initial_seed
00187 
00188     failure=.FALSE.
00189     IF (.NOT.ASSOCIATED(meta_env)) RETURN
00190 
00191     initial_seed = next_rng_seed(error=error)
00192 
00193     DO i_c=1,meta_env%n_colvar
00194        NULLIFY (meta_env%rng(i_c)%stream)
00195     END DO
00196 
00197     ! Each process generates all seeds. The seed generation should be
00198     ! quite fast and in this way a broadcast is avoided.
00199 
00200     ALLOCATE (seed(3,2,meta_env%n_colvar),STAT=stat)
00201     CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00202 
00203     seed(:,:,1) = initial_seed
00204     DO i_c=2,meta_env%n_colvar
00205        seed(:,:,i_c) = next_rng_seed(seed(:,:,i_c-1),error=error)
00206     END DO
00207 
00208     ! Update initial seed
00209     initial_seed = next_rng_seed(seed(:,:,meta_env%n_colvar),error=error)
00210 
00211     ! Create a random number stream (Wiener process) for each particle
00212     DO i_c=1,meta_env%n_colvar
00213        WRITE (UNIT=name,FMT="(A,I8)") "Wiener process for COLVAR",i_c
00214        CALL compress(name)
00215        CALL create_rng_stream(rng_stream=meta_env%rng(i_c)%stream,name=name,distribution_type=GAUSSIAN,&
00216             extended_precision=.TRUE., seed=seed(:,:,i_c),error=error)
00217     END DO
00218     DEALLOCATE (seed,STAT=stat)
00219     CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
00220 
00221   END SUBROUTINE create_wiener_process_cv
00222 
00223 END MODULE wiener_process